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Editorial 


With this first issue of the Society of Petroleum Engineers Journal, the SPE 
reaches another milestone in its history of disseminating technical information to 
engineers of the petroleum industry. Publishing of this new journal represents 
another step taken by SPE and AIME over the years to constantly imptove our 
publications program. 

The: first AIME Transactions volume was published in 1873, but publication 
of periodicals in addition to the bound volume did not start until 1905 with a 
journal called The Bulletin. First appearing on a bi-monthly basis, The Bulletin 
became a monthly publication in 1909. 

In 1919 Mining and Metallurgy replaced The Bulletin and was published through 
1948, serving all members of AIME. It contained no Technical Papers, but included 
many technical articles. During the earlier days of its publication, the Technical 
Papers were issued as separates which members requested by mail. Then in 1934, 
the first of the four Technologies was originated to complement Mining and Metal- 
lurgy and to handle Technical Papers dealing with specific fields of technology. 

Metals Technology was issued first in 1934, Mining Technology in 1937, 
Petroleum Technology in 1938 and Coal Technology in 1946, They were published 
irregularly as bi-monthlies or quarterlies. Then in 1949 publication of Journal of 
Petroleum Technology in Dallas was started and, simultaneously, Journal of 
Metals and Mining Engineering were founded in New York. 

In 1958 the increasing flow of technical material again caused expansion of 
the publications programs. The Society of Petroleum Engineers began to publish 
some papers as separates, with summaries, appearing in Journal of Petroleum 
Technology and the complete text in the annual Transactions volume. At the 
same time, The Metallurgical Society of AIME separated metals Transactions 
papers from Journal of Metals, placing them in a bi-monthly Transactions of the 
Metallurgical Society. Prior to that time, since 1953, Journal of Metals had issued 
three to five supplements per year carrying additional Transactions papers. 

The Society of Mining Engineers decided in 1960 to begin the practice of 
publishing most of their Transactions material in the annual volume only. Now, 
in 1961, the Society of Petroleum Engineers begins publication of Society of 
Petroleum Engineers Journal. 

Looking back over the years in this manner, it is interesting to see how the 
three societies of AIME have developed their publications programs similarly, 
even though they have approached the problem independently, particularly in 
recent years. 

The primary purpose of a professional society is to provide a forum for technical 
information which will aid its members in better fulfilling their jobs, So it is with 
the Society of Petroleum Engineers as we inaugurate this new publication with 
great hopes and expectations for its future. Its success, however, depends not 
upon the ability of the Society to publish the material, but upon the efforts of 
individual members in supplying outstanding Technical Papers for presentation 
here. This is your publication; your contributions will be the foundation for its 
growth and service to the profession. 


R. WILLIAM TAYLOR 
Editor 





SOCIETY OF PETROLEUM ENGINEERS JOURNAL 





Experiments on Mixing During Miscible Displacement 


in Porous Media 


WILLIAM E. BRIGHAM 
JUNIOR MEMBER AIME 
PHILIP W. REED 

JOHN N. DEW 
MEMBERS AIME 


ABSTRACT 


The paper describes experiments on miscible dis- 
placement in various porous media and the results 
of these experiments. Both glass bead packs and 
natural cores were used. Bead diameters varied from 
0.044 to 0.47 mm, and pack lengths varied from 83 
to 678 cm. Natural cores used were Berea and 
Torpedo sandstone. 

By taking samples as small as 0.5 cc and using 
refractive index for analysis, the data on break- 
through curves could be plotted to within +0.5 per 
cent. To plot the data correctly on error function 
paper, a parameter (Vp — V)/\/V was used which 
allowed for the predicted growth of the front as it 
moved past the observer. 

The change in the amount of mixing (length of 
mixed zone) was studied by varying velocity, length 
of travel, bead size, viscosity ratio and pack diam- 
eter. When the displaced material was less viscous 
than the displacing material (favorable viscosity 
ratio), these changes were adequately predicted by 
theory. When natural cores were used, rather than 
glass beads, the amount of mixing was greatly in- 
creased — also qualitatively predicted by theory. 

In experiments with favorable viscosity ratios in 
which the ratio was varied from 0.175 to 0.998, it 
was found that the rate of mixing was changed by a 
factor of 5.7. Thus, the rate of mixing is strongly 
affected by viscosity ratio, even when the theoreti- 
cal error function relationship for mixing is valid. 

Experiments using fluids with viscosity ratios near 
1.0 showed that the instability effects of even a 
slightly unfavorable viscosity ratio (1.002) caused 
disproportionately more elongated breakthrough 
curves than found with a favorable viscosity ratio 
(.998). When the viscosity ratio was as high as 5.71 
these instability effects were much more pronounced, 
as evidenced by the shape of the breakthrough curve. 
The displacements at viscosity ratios above 1.0 no 





Original manuscript received in Society of Petroleum Engi- 
neers Office Nov. 30, 1959. Revised manuscript received Oct. 
13, 1960. Paper presented at AIChHE-SPE Joint Symposium, Dec. 
6-9, 1959, in San Francisco. 








MARCH, 1961 





CONTINENTAL OIL CO. 
PONCA CITY, OKLA. 


longer followed the theoretical error function curve. 

It was found that at reservoir rates in natural res- 
ervoir materials the rate of mixing (as measured by 
the dispersion coefficient K) was considerably 
higher than in glass bead packs. The higher mixing 
rate is caused by the inhomogeneity of reservoir 
rock as compared to glass beads. The amount of 
inhomogeneity is expressed in terms of a mixing 
coefficient a. 

Data on length of mixed zone vs velocity showed 
that there exists a velocity at which the zone length 
is a minimum. At this velocity, diffusional mixing 
contributes only a small fraction of the total mixing. 


INTRODUCTION 


Data have been taken on the amount of mixing 
between two miscible fluids during the displacement 
of one fluid by another using various systems of 
porous media and various fluids. The porous media 
were packs of relatively uniform spherical glass 
beads (which form an ideal medium for flow studies) 
and consolidated sandstone cores. The effects of 
displacing rate, bead size, length of travel, diameter 
of bead pack and viscosity ratio were investigated. 


EQUIPMENT AND EXPERIMENTAL PROCEDURE 


The glass beads were packed in Lucite tubes 
which had an inside diameter of 3.19 cm and a wall 
thickness of 0.64 cm, except for those packs used 
to study diameter effects. Porous stainless steel 
discs were set in the ends of the tube to contain 
the beads. Minnesota Mining and Manufacturing Co. 
*‘Superbrite’’ beads were used. The consolidated 
cores were mounted in Lucite tubes, which had been 
heated to the softening point and bonded to the cores 
using external pressure. Data on the bead packs 
and cores are given in Table 1. The prefix numbers 
on the bead packs (107, 113 and 117) are the Min- 
nesota Mining and Manufacturing code numbers for 
the particular bead size used. 

Fluid movement at any desired rate was obtained 
by a Zenith gear-type pump driven through a variable 
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speed transmission. The bead packs (or cores) were 
mounted vertically with the displacing oil injected 
at the bottom. Oils of matched density were used. 
Oil compositions are given in Table 1. The diffusion 
coefficient of these fluids was found experimentally 
to be 3.2 x 10-6 cm2/sec. The proportion of each 
oil in the effluent was determined from refractive 
index (R.I.) measurements since one oil was mainly 
Soltrol (an isoparaffinic oil with a low R.I.) and 
the other was mainly commercial kerosene with a 
high R.I. For convenience in identification and 
visual observation of the fluid displacement, one 
of the pair of oils was dyed red with a dye which 
did not change the refractive index. 

In preparation for a run, a pack was saturated 
with oil to be displaced. The oil-saturated pack was 
inverted, the header filled with displacing oil and 
the tubing carrying displacing oil was connected. 
The pack was re-inverted so that the inlet was at 
the bottom, the outlet header was emptied and the 
run started immediately. In this way, an exact ma- 
terial balance could be made on the effluent. The 
effluent was caught in a graduate until the approx- 
imate arrival of the mixed zone. Incremental samples 
of suitable size were then taken as the effluent 
accumulated in the header. Samples were taken with 
a hypodermic syringe so that each sample emptied 
the header; thus, each sample was exactly repre- 
sentative -of the oil passing the end of the pack 
during the small time interval of that sample. The 
refractive index was measured and the proportions 
of displacing and displaced oil determined from an 
experimental calibration curve. By taking samples 
as small as 0.50 cc, the data on breakthrough 
curves could be plotted within +0.5 per cent. 


EXPERIMENTAL RESULTS 


EFFECT OF LENGTH OF TRAVEL 
Taylor’s! theory of displacement in capillary 





lReferences given at end of paper. 


tubes, Scheidegger’s? statistical theory of porous 
media, Keulemans’? ‘teddy diffusion’? theory and 
Frankel’s4 ‘'stagnant pockets’’ theory all predict 
that longitudinal dispersion will be governed by 
the following equation. 

oC . K 02C 

ot Ox 42 

Eq. 1 is identical in form to the well known dif- 
fusion equation; however, the rate of dispersion or 
mixing of the flood front is governed by a disper- 
sion coefficient (K) rather than the Fick diffusion 
coefficient (D), as in the diffusion equation. Also, 
since the flood front is moving through the porous 
medium, the space variable used is the distance 
from the midpoint of the flood front (x, ) rather than 
the distance from the inlet end of the porous medium 
(x). The variables x and x, are simply related as 
follows: 


ee ae eee ore ae ee 


= ass L 
Ry RE EE ae ak eee xe | 
The solution of Eq. 1 with appropriate boundary 
conditions is well known in the following form. 


= ] — erf 
y2 \1 art (3*te) hg - 


where 
2 
erf (é) ~ Vir Fefae : 


Eq. 3 defines the well known ‘‘S’’ shaped curve 
of the error function integral, and a plot of this 
equation is a straight line on arithmetic probability 
co-ordinate paper. 

The argument of the error function in Eq. 3 
(x,/2\/Rt) indicates that at a constant rate of flow 
and with a constant dispersion coefficient the 
spread of the mixed zone will be proportional to 
the square root of the distance traveled. It is of 
interest to determine the validity of this ‘‘square 
root law’’; however, first it is necessary to change 








TABLE 1 — POROUS MEDIA AND FLUIDS 

















Pack or Length Diameter Pore Vol. Porosity: Permeability Avg. Bead Dia. 
Core No. Lem) _ fem) _ ae. (per cent) _(darcies) _ (mm) 
117-1 83.2 3.19 221 33.2 0.65 044 
113-1 83.3 3.19 230 34.6 7.1 - 100 
107-1 83.5 3.19 215 32.2 115 -470 
113-2 167.2 3.20 500 37.1 7 - 100 
113-3 83.4 14.00 4615 36.0 Y 4 -120 
113-4 83.5 1.30 47 37.0 7 -100 
113-6 678 3.19 1946 35.8 7 -100 
B-1 82.2 7.60 725 20.0 0.300 Berea Sandstone 
-042 (Calculated) 
T-1 82.2 4.88 346 22.5 0.258 Torpedo Sandstone 
-036 (Calculated)’ 
Oils for Favorable Viscosity Ratio Runs: 
Displaced Oi! — 89.4 vol. per cent Soltrol, 10.6 per cent Carbon Tetrachloride 
Density 0.840 gm/cc, Viscosity 1.25 cp at 25°C, 
Displacing Oil — 50 vol. per cent Mineral Oil, 50 vol. per cent Kerosene 
Density 0.840 gm/cc, Viscosity 7.14 cp at 25°C. 
Oils for Runs at Approximately Equal Viscosity Ratios: 
Clear Oil — 94.9 vol. per cent Soltrol, 5.1 per cent Carbon Jetrachloride 


Density 0.795 gm/cc, Viscosity 1.296 cp at 25°C, 


— 94.8 vol. per cent Kerosene, 5.2 per cent Heptane 


Density 0.795 gm/cc, Viscosity 1.299 cp at 25°C, 


Red Oil 
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the equation variables somewhat. Eq. 3 defines the 
concentration profile in the porous medium at a 
given time, while the experimental data give the 
concentration at a given point as the constantly 
growing mixed zone moves past the observer. The 
following substitutions can be made: 


Vp-V 
x,=L _ a a ae ae ee ee er ee ee (4) 
Dp 


(5) 


HT] 
=< 
~~ |< 
Sensei 
° 
Ln | 
It 
|x 


Thus, the argument of the error function becomes, 





x L(V,-V) E 
es oe oe  (, .  . e 
2VKt 2VTKV,V 2VKTV, 
where 
V,—-V 
v=. 


VV 

The error function parameter U serves two pur- 
poses. First, it accounts for the predicted growth 
of the mixed zone as it passes an observer sta- 
tioned at a fixed point. Second, with displacements 
varying only in length of travel, plots of U vs con- 
centration should fall on the same line if the ‘‘square 
root law’’ holds. 

As atest of the square root law, runs were made 
at approximately the same rate (5.2 to 6.3 x 1073 
cm/sec) at a viscosity ratio of 0.175. Three bead 
packs of 3.19-cm ID were used — Pack Nos. 113-1, 
2 and 6. The lengths of these packs were 83.3, 
167.2 and 678 cm, respectively. The results of 
these runs are plotted in Fig. 1. As can be seen, 
the data fall on a straight line in the concentration 
range of 5 to 80 per cent displacing fluid. This 
shows that the theoretical equation for the com- 
position profile (the error function) is valid over 
most of the concentration range, although at high 
concentrations the profile is somewhat longer than 
predicted. Further, since all the data fall on the 
same straight line, the square root law is also 
valid. 

It appears appropriate at this point to define the 
relationship between the error function parameter U 
and the rate of growth of the mixed zone which is 


© 63.3 cm BEAD PACK, NO 113-1 
& 167.2 cm BEAD PACK NO. 113-2 
OG 678 cm BEAD PACK NO. 113-6 





% DISPLACING FLUID 


FIG. 1— EFFECT OF LENGTH OF TRAVEL ON ERROR 
FUNCTION PLOT. 





MARCH, 1961 








characterized by the dispersion coefficient K. 
Taylor! has shown that K can be related to the length 
of the transition zone between any two specified com- 
positions. For instance, if xjq is defined as the 
distance to the point of 10 percent displacing fluid 
and xgq is the distance to the point of 90 per cent 
displacing fluid, the length of transition zone so 
defined is x9q— x19, and the dispersion coefficient 
is defined as follows. 


2 
oe Eea eee a ee ee eee 
t 3.625 


However, for our purposes it is necessary to relate 
the dispersion coefficient K to the error function 
parameter U. The following method is used. The 
best straight line is drawn through the data points. 
Values of U are read from this line at the 10 per 
cent (U,9) and 90 per cent (Ugg) concentration val- 
ues. The dispersion coefficient is then defined as 
follows. 


2 
a L (Ugo ~ U10) 
Vp'T 3.625 





. (8) 


Eq. 8 is identical to Eq. 7, as can be shown by 
substituting Eqs. 5 and 6 into Eq. 7. If a person 
prefers to use other concentration values, such as 
20 and 80 per cent, the form of the equations would 
be the same as Eqs. 7 and 8, but a different con- 
stant would be necessary (e.g., 2.380 for 20 and 
80 per cent concentrations), These constants can 
be obtained from any standard table of error inte- 
grals. 


EFFECT OF VISCOSITY RATIO 

Data were run in the 113-1 bead pack at viscosity 
ratios of 0.175, 0.998, 1.002 and 5.71 at a constant 
velocity (= 6.0 x 1073 cm/sec). The effluent com- 
position curves are shown in Figs. 2 and 2A. In 
Fig. 2 only the 5.71 and 0.998 viscosity ratio data 
are shown due to crowding of the other curves near 
1.0 pore volume injected. In Fig. 2A the central 
portion of the effluent composition plot of Fig. 2 
has been greatly expanded, and the data at 0.175 
and 1.002 viscosity ratios have also been included. 
Of course, on this expanded plot only a small por- 
tion of the 5.71 viscosity ratio curve can be seen. 

Fig. 2 best shows the marked difference when 
the viscosity ratio of 0.998 (a favorable viscosity 
ratio) is compared to a ratio of 5.71 (an unfavorable 
viscosity ratio). The favorable ratio curve is the 
typical ‘‘S’’ shaped curve of the error function, as 
defined by Eq. 3. The unfavorable viscosity curve 
is greatly extended, no longer follows the theoreti- 
cal ‘‘S’’ shaped curve and exhibits a number of 
bumps and irregularities in composition. The irreg- 
ularities are caused when ‘‘fingers’’ of displacing 
fluid run ahead of the bulk of the flood front. These 
fingers could be easily observed since the displac- 
ing fluid was dyed red; also, they have been reported 
visually by Blackwell, et al.5 

In Fig. 3, the data of Fig. 2A have been plotted 
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on probability co-ordinates. The most noticeable 
phenomenon is that increasing the viscosity ratio 
also increases the value of Ugg — Ujo and, thus, 
increases K and the rate of dispersion. For instance, 
compare the data at viscosity ratios of 0.175 and 
0.998. Both of these are favorable ratios and both 
follow the theoretical error function of Eq. 3. But, 
it is evident that the change in viscosity ratio 
caused an increase in mixing rate (shown through 
Ugg — Uyo). This phenomenon has not been ade- 
quately studied to date, so no generalization can be 
made about the effect of viscosity ratio on mixing 
rate. But this points out the need to view with care 
data from different authors who may be using dif- 
ferent viscosity ratios. 

Another important effect of viscosity ratio is found 
by comparing the 0.998 data with the 1.002 data. 
The instability effects of only a slightly unfavorable 
viscosity ratio (1.002) cause a disproportionately 
large change in the effluent concentration curve. 
Also, along with a greater amount of mixing, the 
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PORE VOLUMES INJECTED 


FIG, 2— EFFECT OF VISCOSITY RATIO ON EFFLUENT 


CONCENTRATION CURVE, BEAD PACK NO. 
113-1. 


% DISPLACING FLUID 


PORE VOLUMES INJECTED 


FIG. 2A—EFFECT OF VISCOSITY RATIO ON EFFLU- 
ENT CONCENTRATION CURVE, BEAD PACK 


NO. 113-1. 











unfavorable viscosity ratio data no longer plot a 
straight line on the probability co-ordinates of Fig. 
3. The 5.71 viscosity ratio data are not plotted on 
Fig. 3, but these data exhibit even more curvature. 
So Eq. 3 cannot be considered valid for displace- 
ments at a viscosity ratio greater than unity. 


EFFECT OF PACK DIAMETER 

Displacements were run in 14.00-cm (113-3), 3.19- 
cm (113-1) and 1.30-cm (113-4) packs to determine 
the effect of over-all pack diameter. The beads 
used in all three packs were 0.100-mm average 
diameter, and the viscosity ratio used was 0.175. 
As can be seen in Fig. 4, the 3.19 and 14.00-cm 
packs show almost identical breakthrough curves; 
but, the 1.30-cm pack definitely exhibits more dis- 
persion. There are two possible reasons for the’ 
greater mixing found in the small-diameter pack. 
One reason can be deduced from the nature of the 
mixing theories in porous media.2,3,4 All these 
theories assume a cross section large enough that 
boundary or wall effects are negligible, and it is 
likely that this assumption is not valid in the 1.30- 
cm pack. Another reason for greater dispersion may 





% DISPLACING FLUID 


FIG. 3—EFFECT OF VISCOSITY RATIO ON ERROR 
FUNCTION PLOT, BEAD PACK NO, 113-1. 
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= 3.19 CM DIA. BEAD PACK — 113-1 


* 13.CM DIA. BEAD 
PACK - 113-4 & = 14.00 CM DIA. BEAD PACK-113-3 





1.01 1.02 1.03 1.04 1.05 1.06 


0.94 095 036 100 ° 
PORE VOLUME INJECTED 


FIG. 4—EFFECT OF BEAD PACK DIAMETER ON EF- 
FLUENT COMPOSITION CURVE. 
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be the difficulty in achieving uniform packing in a 
small tubing; but, whatever the reason, the data of 
Fig. 4 emphasize the need for viewing with caution 
any data acquired from small-diameter packs (less 
than approximately 3-cm diameter). 


EFFECT OF FLOW VELOCITY 

Keulemans3? has suggested that one source of 
dispersion at higher flow rates is due to the separ- 
ation and rejoining of channels of flow in the porous 
medium and has termed this ‘‘eddy diffusion’. If 
one considers a particular porous medium, the amount 
of dispersion introduced with each centimeter of 
distance traveled is the same, regardless of the 
velocity of flow through that medium. The dispersion 
is-entirely due to the amount of separation and re- 
joining that is inherent in the medium itself. How- 
eve:, the dispersion coefficient (K) is a measure of 
the rate of dispersion, not the amount, and thus 
would be proportional to velocity. Also, when con- 
sidering two similar porous media differing only in 
size of flow channels and particles, the amount of 
dispersion is directly proportional to particle size. 
So, for Keulemans’ ‘‘eddy diffusion’’, 


CEA 44.444 44 64 # eo 6 oe @ 
where A is a term which expresses the amount of 
irregularity of packing and, thus, the amount of 
eddy diffusion characteristic of the particular porous 
medium. It is often desirable to express variables 
in dimensionless form, and Eq. 9 suggests the fol- 
lowing form. 

Km 

Dp 6» 

Taylor’s capillary tube theory and Frankel’s 
‘stagnant pockets’’ theory both reduce to the fol- 
lowing proportionality, 


(10) 


Tu)2 
ol (11) 
which, again in dimensionless form, becomes 
K (ru 
D ~\pD (12) 


By referring to Expressions 10 and 12, it would 
appear reasonable to plot 7u/D vs K/D, both for the 
purpose of correlation of data and to test the valid- 
ity of these theories. Accordingly, a plot of all the 
data has been made in Fig. 5 using K/D and 7u/D 
as co-ordinates, as done by Blackwell, et al.5 A 
log-log co-ordinate system is preferable here be- 
cause of the wide range of variables covered and 
because 7u/D is expressed as a power function in 
the above propo? ‘onalities. 


EFFECT OF FLOW VELOCITY — LOW FLOW RATES 
To understand the amount of mixing which occurs 
at low velocity, one should consider the limiting 
case where velocity is zero. At non-flow conditions, 
the only case of mixing should be ordinary molecular 
(Fick) diffusion, and the dispersion coefficient K 
will be constant. In Fig. 5, the data on 0.100-mm 
beads approaches this condition at a value of 7u/D 
of about 0.2. Thus, at low rates the equation for the 
dispersion coefficient in glass bead packs becomes 


MARCH, 1961 








K 
a ee ee (13) 
Blackwell, et al> and Carman® have also found 
this relationship at low flow rates for spherical 
glass bead packs. The number 0.7 should not be 
considered a universal constant, for diffusion in a 
porous medium is directly analogous to the electri- 
cal resistivity of the medium. Thus, K/D is related 
to the formation factor F and porosity ¢, and a more 
exact equation at low rates is 
K 1 
The term 1/Fq@ will commonly vary between 0.15 
and 0.7, depending on the lithology of the porous 
medium. 7 


EFFECT OF FLOW VELOCITY — HIGH FLOW RATES 

At high flow rates, as stated, Keulemans’ theory 
predicts that the term K/D is proportional to 7u/D 
to the first power. The other theories predict a 
second power dependence. The data of Fig. 5 show 
that none of these theories is exact but, rather, that 
the result lies somewhere between them. The expo- 
nent was found experimentally to be 1.24, 1.20 and 
1.19 in the 0.470-mm beads, 0.100-mm beads and 
the Berea core, respectively. These data compare 
closely with the 1.17 power reported by Blackwell, 
et al,> and the 1.20 power shown by Aronofsky and 
Heller8 on data reported by von Rosenburg.9 

From these data it appears that, for glass bead 
packs and fairly homogeneous sandstone, a reason- 
able value to use for the exponent is 1.20. Thus, we 
get the empirical equation 


1.20 
5 ~a(F) Dk het eeu 


The exponent (1.20) may not be a constant, valid 
for all porous media, since it has only been tested 
for bead packs and sandstone. It lies between the 
1.00 and 2.00 limits of the mixing theories, so it 


BEADS .470 MM DIA. U/ 2 * 0.175 
BEADS . 100 MM DIA. 4 /Mg © O75 


BEADS . 044 MM DIA. 4, /ug = 0.175 
BEADS 1 100 MM DIA. U,/uz = 0.998 
BEREA .300 MO PERM. U;/lp= 0.175 
TORPEDO 258 MD PERM. 4h /lg=0.175 


470 MM DIA 


iu 
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FIG. 5—EFFECT OF RATE AND TYPE OF POROUS 
MEDIUM ON DISPERSION RATE. 











appears that both eddy diffusion of Keulemans’ 
model and statistical dispersion of either the 
Scheidegger or Frankel type play a part in deter- 
mining the rate of mixing. Also, since the result is 
close to the first power, eddy diffusion must be the 
dominant factor in these media. 

Eq. 14 for low rates can be combined with Eq. 15 
for high rates to give a single equation expressing 
the dispersion coefficient as a function of velocity 
and lithology. Thus, we get 


ee ae ecules oe Jeena 


It will be noted that no attempt has been made to 
give an exact numerical value for the mixing coeffi- 
cient (a) in these equations. The reason is that ais 
not a constant, as implied in other work® but, rather, 
varies widely depending on the geometry of the 
porous medium and the viscosity ratio of the flow- 
ing fluids. A further discussion on the effect of 
porous medium type on a is found in the next sec- 
tion, and a comparison of a values has been made 
in Table 2. The authors feel that, with more exper- 
imental evidence, the mixing coefficient (a) could 
be used with success to characterize porous media, 
for ais a measure of the inhomogeneity of a medium 
in a dynamic situation, while the present pore size 
distribution techniques are based on non- flow 
measurements. 

Upon studying Eq. 16, one observes an interest- 
ing relationship between the dispersion coefficient 
K and the molecular diffusion coefficient D. At 
high flow rates with all other variables held con- 
stant, K is inversely proportional to D®?; this 
means that, if the fluids had a larger diffusion co- 
efficient, the amount of mixing would be slightly 
less. No attempt has been made to study this pos- 
sibility here, but the data of Handy! at least 
qualitatively substantiate that mixing rate is de- 
creased somewhat when using a fluid with a greater 
diffusion coefficient. 


EFFECT OF POROUS-MEDIUM TYPE 
Keulemans’ eddy diffusion term includes the co- 
efficient } which expressed the inhomogeneity of 
the flow channels. Let us assume we have two 
porous media which are geometrically similar with 
different particle sizes, such as two identically 
packed bead packs differing only in diameter of 
beads. These packs would have identical values 
for A according to Keulemans’ theory. Scheidegger’s 
theory contains a dispersivity term a which from 
its definition would also be constant under the 
afore-mentioned conditions, and Frankel’s stagnant 





TABLE 2 — EFFECT OF VISCOSITY RATIO AND POROUS 
MEDIUM ON @, FOR K/D = a (tu/D) 1-20 








Porous Medium Viscosity Ratio a 
-044-mm Beads ATS 0.69 
-100-mm Beads 0175 0.49 
- 100-mm Beads -998 2.78 
-470-mm Beads 0175 0.30 
Torpedo Sandstone 175 23.2 
Berea Sandstone 175 53 


6 


pockets theory gives the same result. However, in 
practice, A (and, thus, a and a) has experimentally 
been found to increase with finer particles. Klinken- 
berg and Sjenitzer 11 report values of A increasing 
from one to eight as particle diameters drop from 
1.0 to .05 mm. The data of Fig. 5 and Table 2 on 
0.470, 0.100 and .044-mm beads also point out this 
trend, for the values of a increase as the bead 
diameters decrease. 


Why then is a not constant? Keulemans suggested 
that this increase with small particles is due parti- 
ally to the difficulty in packing smaller particles 
uniformly and partially to the greater variation in 
particle size usually found with finer particles. 
That is, the larger a is caused by more inhomoge- 
neity in the porous medium. The data reported here 
from glass bead packs corroborates this theory, for 
the small beads have a greater variation in particle 
size than the large beads. 

Also, we can compare the values for a reported 
here to those cited by Blackwell5 and Klinkenberg. 11 
The a values found in this paper are considerably 
smaller than found by these two authors. Part of 
this difference is caused by viscosity ratio (ex- 
plained earlier), but even the 0.998 viscosity ratio 
data give less mixing than reported by Blackwell 
and Klinkenberg. Again, this is caused by the nar- 
row particle-size range of these beads — that is, 
by the more homogeneous porous media. 

To plot the Berea and Torpedo sand data on Fig. 
5 for comparison with the bead pack data, an as- 
sumption must be made on average particle radius 
to be used for these natural porous media. The fol- 
lowing equations were used to calculate average 
pore radii and particle radii. 


zy ~ /8k r2 a4 
Cg OP gk ge ee ee 


7 = 3r, eee ee ee a ee ee ee 
where 
7 = 2.0 


Using these equations, the average sand-grain 
radii calculated 2.0 x 10-3 cm for the Berea sand- 
stone and 1.8 x 1073 cm for the Torpedo sandstone. 
These average grain radii are, of course, not exact, 
but they are close enough to show without doubt 
the tremendously increased importance of the a 
term (the inhomogeneities) in naturally occurring 
porous media when the data are compared in Fig. 5 
to the homogeneous bead pack. Handy’s!© data 
obtained from a Boise consolidated sandstone show 
the same results. In Table 2, the a values for the 
bead pack data and the natural core data may be 
easily compared. 

To more readily see the effect of inhomogeneities 
on rate of mixing at reservoir rates, let us consider 
a displacement at 3.0 x 10-4 cm/sec (0.85 ft/day). 
In the .100-mm glass bead pack, the dispersion rate 
was controlled almost entirely by molecular diffu- 
sion. But, in the Berea core used in this study, the 
dispersion rate was 9,2 times as great as in the 
glass bead pack because natural sandstone is more 
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heterogeneous. This result is reflected directly in 
the values for a found in Table 2. 


EFFECT OF VELOCITY ON ZONE LENGTH 
It is of considerable interest to know the effect 
of velocity on the length of the mixed zone. Con- 


sider Eqs. 7 and 16. We can substitute 5¥ for t in 
p 


Eq. 7, and then solve these equations simultaneously 
for x99 — x39 as a function of velocity. The result 


is 
(‘20-210 a 7, BLY «(22° (BLY) 
3.625 Fe uv, D uV, 
+ ae SS ee ee (19) 
A glance at Eq. 19 will show that there exists a 
velocity at which the mixed zone length is a mini- 
mum. Eq. 19 can be differentiated with respect to 
velocity and set equal to zero to solve for this 
velocity. 


iia 1 1/1.20 
-Ftae) ate ene San, 


Using the data for the 113-1 bead pack (0.70 for 
1/F¢ and 0.49 for a), the velocity for minimum mix- 
ing is found from Eq. 20 to be 3.3 x 10-3 cm/sec. 
Fig. 6 shows the data on velocity and zone length 
for this pack along with the calculated minimum. 

Another method of fixing this ‘‘minimum’’ velocity 
is to draw a 45° line on Fig. 5 tangent to the curve. 
The point of tangency corresponds to the minimum 
point. By this procedure it is apparent that this 
velocity will be lower as the porous medium be- 
comes more heterogeneous. For instance, assuming 
1/F@ equals 0.50 for the Berea sandstone and a 
equals 53, the velocity for minimum mixing in Berea 
is calculated to be 1.2 x 10-4 cm/sec. Also, it can 
be seen from the K/D value that only about 20 per 
cent of the total dispersion coefficient is caused 
by diffusion at the ‘‘minimum mixing’’ velocity. 





SUMMARY AND CONCLUSIONS 


This investigation has found that the viscosity 
ratio during displacement greatly affects dispersion 
rate in the following ways. 


LENGTH OF MIXED ZONE (X90-Xi0)(cm) 





VELOCITY, u ( cm/sec) 


FIG. 6—EFFECT OF VELOCITY ON ZONE LENGTH, 
BEAD PACK NO. 113-1, VISCOSITY RATIO 
0.175. 
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1, When a favorable viscosity ratio was changed 
by a factor of 5.7 (0.175 to 0.998), the dispersion 
coefficient was also changed by a factor of 5.7. 
Although this seems to suggest a linear relation- 
ship, there are not enough data available to sub- 
stantiate such a conclusion. This should be in- 
vestigated further. 

2. When changing only slightly from a favorable 
(0.998) to an unfavorable (1.002) viscosity ratio, 
the mixing rate was disproportionately increased 
and the predicted error function curve was no longer 
valid. 

3. When the viscosity ratio was 5.71, the effluent 
concentration curve was greatly extended, and 
there was both visual and analytical evidence of 
distinct ‘‘fingers’? of displacing fluid running 
ahead of the average displacement front. 

4, At viscosity ratios above 1.0, the theoretical 
error function curve no longer was valid. 

The following conclusions can be drawn for dis- 
placements at a favorable viscosity ratio (viscosity 
ratio < 1.0). 

1. The ‘‘square root law’’, which concludes that 
the amount of mixing is proportional to the square 
root of the distance traveled, is valid. 

2. Diameter of the pack can be an important 
variable, and data from small-diameter packs should 
be treated with caution. Data obtained using .100- 
mm diameter beads were consistent for packs 3.2 
cm in diameter and larger. 

3. At very low rates of flow, the dispersion co- 
efficient K is directly proportional to the Fick 
diffusion coefficient, as expected. 

K 1 


D Fd¢* 

4. At high rates of flow, the dispersion coeffi- 
cient is characterized by the following relation for 
sandstones and glass bead packs. 


K ais Tu 1.20 

D D 

5. a in this equation is not a constant but, rather, 
is a function of the viscosity ratio of the flowing 
fluids and the inhomogeneity of the porous medium. 

6. In natural sandstones at reservoir rates, the 
dispersion coefficient K is proportional to (u)!+2°, 
as in the above equation, and the rate of dispersion 
is considerably higher than found in glass bead 
packs. 

7. The length of the mixed zone is a function of 
velocity. The zone length increases at very low flow 
rates and very high flow rates, and there exists a 
velocity at which the zone length is a minimum. At 
this velocity, diffusion contributes only a small 
fraction of the total dispersion coefficient. ' 


NOMENCL ATURE 


= concentration (cc/cc) 

= time (seconds) 

x 1 = distance from midpoint of flood front 
(cm) 

dispersion coefficient (cm2/sec) 


mm O 
mY 


za 
Il 






























= Fick diffusion coefficient (cm*/sec) 
x = distance from inlet end of porous me- 
dium (cm) 
u = average pore velocity (cm/sec) u=L/T 
length of porous medium (cm) 
T = time required to inject or produce one 
pore volume of fluid (seconds) 
Vp = pore volume of porous medium (cc) 
V = volume of fluid recovered at time of 
sample (cc) 


| ad 
i} 


U = error function parameter(V, — V/V 

7 = average radius of glass beads or sand 
grains (cm) 

E = Keulemans’ eddy diffusion coefficient 
(cm2/sec) 

A = Keulemans’ packing coefficient 

F = formation resistivity factor 

¢ = fractional porosity 


mixing coefficient 
Scheidegger’s dispersivity 
= average pore radius (cm) 

= permeability (darcies) 


7, 


ws aR 
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Viscosity 
Ratio = the ratio of the displaced phase 
viscosity to the displacing phase 
viscosity 
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Stability Theory and Its Use to Optimize 


Solvent Recovery of Oil 


RICHARD L. PERRINE 
MEMBER AIME 


ABSTRACT 


This paper shows how stability theory can be 
used to optimize solvent recovery of oil. Applica- 
tion of the theory leads to definition of the limiting 
conditions required for stable displacement to 
occur. One of these conditions is that the size of 
a solvent ‘‘slug’’ must at least equal a prescribed 
minimum. The criteria to be satisfied are miscibil- 
ity and stability. Stability implies that no viscous 
fingering will occur and that mixing will be caused 
only by a dispersion process. Miscibility implies 
that complete recovery will be obtained from the 
swept regions. Thus, for any specified set of res- 
ervoir conditions, an optimum use of solvent is 


defined. 


INTRODUCTION 


The early outlook for solvent flooding as a means 
to increase cil recovery was very favorable. Some 
laboratory results indicated that a ‘‘slug’’ con- 
taining perhaps 2 to 3 per cent of a hydrocgrbon 
pore volume could be successful. However, pther 
data suggested as much as 30 per cent would be 
required. The difference is of considerable economic 
importance. A theoretical explanation for these 
divergent results has been advanced by the author.!»2 
Stability theory defines conditions for two distinct 
flow regimes. In stable displacement, solvent and 
oil become mixed by a dispersion process, and the 
solvent requirement is small. On the other hand, 
unstable displacement can degenerate into viscous 
fingering. The practical result is a considerable 
increase in the extent of mixing and, thus, in the 
solvent required. 

The present paper shows how to use stability 
theory to optimize solvent recovery of oil. Oil is 
usually considered to be displaced by a small amount, 
or a slug, of solvent. Gas in turn follows solvent, 
Any fingering will permit contact of nearly solvent- 
free gas and oil, leading to immiscibility and re- 
duced recovery. Thus, our optimum is defined by 
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two restrictions — stable flow and miscibility. 
Although an immiscible solvent process could prove 
more profitable than any miscible process, such 
processes lie outside the selected definition of an 
optimum. The usefulness of computations based on 
the present work will depend on the validity of the 
general theory. Because the impor-ant idealizations 
required tend to minimize the predicted solvent re- 
quirement, the results should state correctly the 
limiting conditions required. Therefore, theoretical 
calculations should have practical value. 

This paper first reviews several principles that 
result from stability theory. The following section 
describes use of these ideas in relatively simple 
systems which symmetry makes appear one-dimen- 
sional. Of particular importance is the step-by-step 
design of a stable minimum-solvent slug. The use 
of stability theory in cases of multi-dimensional 
displacement is discussed and, finally, suggestions 
are given for simplified practical application of the 
theoretical results. 


PRINCIPLES OF STABILITY THEORY 


The use of perturbation methods to derive stabil- 
ity conditions has been shown by the author.! The 
discussion is presented in terms of the fractional 
concentration of solvent in the single hydrocarbon 
phase. Quite generally, any disturbance can be 
represented as a Fourier series, of which only the 
least stable term need be considered. Denoting 
this term at time ¢t =0 by cy, 


c, = ¢, cos axL~! cos fy! cos nmL“"....(1) 


n n 
The initial growth of this least stable term is cal- 
culated according to 


Ap 0H FR) ce nnnnmnncninnmnienmmaninmnl 8) 


where 
L-?|a2p, + (6? + n?n?)D | 


i 


y 


+ 


-2 oc 2 C7 
U9 [anne - (a + BS] 6, 


+ 
4 
oO 
m 
be 
Q 
= 
| 
t 
Res) 
i) 
+ 
= 
i] 
> 
Lo] 
4 
a | 
aR 
_ 


Oo Oe ee Oe ee es eee ees OSES SESS OSE EE SESE EE SS ESSE ESSE OS OE SESE OS OSES OSES ESOS SEES eS 








and the functions are defined by 


_dinyw,k dp g sin 6 
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g,=2 dp g cos @ 
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The least stable term has the highest growth rate, 
yet stability requires that no growth occur. This 
means that the stability coefficient y must be zero 
or positive, and Eq. 3 shows how this requirement 
restricts displacement conditions. 

The derivation assumed flow in the x-direction 
within a reservoir of large x and y dimensions and 
limited thickness L in the z-direction. Dy is the 
axial and D3 the lateral dispersion coefficient.>-> 
The constant, steble volumetric flow rate per unit 
of cross-sectional area and per unit of porosity is 
ug. Flow symmetry maintains solvent concentration 
c for stabie flow independent of y. Viscosity p is 
for a fluid mixture of solvent concentration c. From 
Eq. 1, a, 8 and n are Fourier expansion parameters 
depending on the initial disturbance. Other symbols 
are k for permeability, ¢ for porosity, p for density, 
g for the acceleration of gravity and @ for the dip 
angle of the formation. For downward tilt in the di- 
rection of flow, 6 is negative. 

Consider what happens as solvent moves through 
a reservoir displacing oil. Dispersion will cause 
mixing and a gradual transition from oil to solvent. 
Values of y will vary within the solvent bank de- 
pending on the relative importance of gravity, vis- 
cosity and dispersion. So long as any part of the 
system is unstable, viscous fingers can grow to 
dominate flood behavior. The problem is to keep 
the most negative value of y zero, forcing all others 
to be positive. No real system is simple. However, 
several useful generalizations can be derived and 
are presented in the following discussions. 


THE PRINCIPLE OF GRAVITY STABILIZATION 

The most important stabilizing mechanism is 
gravity segregation. The function g, in Eq. 4 is 
positive when gravity is dominant, but it becomes 
negative when gravity is subordinate to viscous 
forces. Thus, it measures the relative importance 
of gravity. This assumes ‘‘normal’’ solvent-oil 
density, viscosity and displacement conditions. 
Gravity stabilization can be achieved by control- 
ling rate. 


THE MINIMUM-SLUG-SIZE PRINCIPLE 

The stability limit is marked by y =0. This same 
condition defines the minimum solvent slug for 
stable displacement. Take the simple example in 
which concentration is independent of z and make 
the assumption that a is negligible. Inserting these 
conditions together with y = 0 and re-arranging 


terms in Eq. 3, ac ( p2 + n*7?)D, 





(5) 
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The asterisk indicates a critical, or limiting, re- 
lationship. Stability results whenever either one 
of the following occurs. 

1. gy is positive and (a) the concentration gra- 
dient has any negative value or (b) the concentration 
gradient has a positive value no larger than that 
indicated by Eq. 5. 

2. g1 is negative and (a) the concentration gra- 
dient has any positive value or (b) the concentration 
gradient has a negative value no larger in magni- 
tude (no more negative) than indicated by Eq. 5. 

While a solvent slug should be kept small, at high 
rates stability and an abrupt change from oil to 
solvent may be incompatible.© Maximum gradients 
for a given rate will depend on system properties. 
SOURCES OF REQUIRED DATA 

Most information required to apply Eqs. 3 through 
5 is readily obtained. Although the dispersion co- 
efficients Dj and D3 are new parameters, measure- 
ment techniques have been developed.?+7-9 These 
coefficients characterize the mixing capability of 
a porous medium in much the same way that per- 
meability denotes flow capacity. Because at very 
low rates the lateral «ispersion process will reduce 
to molecular diffusion, these coefficients must be 
measured under appropriate experimental conditions. 

The remaining new parameters arise from the 
Fourier series representation of possible disturb- 
ances in solvent concentration. They are denoted 
by the letters a, B andn. As shown by Eq. 1, these 
parameters reflect the scale, or ‘‘wave length’’, 
of disturbances propagated as a result of inhomo- 
geneity in reservoir properties. The principal cause 
would be permeability inhomogeneity. Means for 
proper determination of a, 8 and n will require 
research into the variation characteristic of porous 
media. Such studies have yet to be made. 

From theory,! however, together with even cur- 
sory observation of natural oil-reservoir rocks, we 
can draw several useful conclusions. The first is 
that very frequently a, and perhaps also f, can be 
neglected. The most important quantity is n. It 
must be an integer and must be greater than or 
equal to one. Thus, conservative practice would 
dictate the values a = 0, 8B = 0 and n = 1 when 
estimates must be used. 

Before assuming that this crude approach will 
always be necessary, consider possible sources of 
information. What we require is a measure of the 
greatest distance between successive extreme 
values of permeability characteristic of the forma- 
tion. In the vertical direction, for example, suppose 
that permeability maxima occur (at most) each L’ 
units apart. Then, it is reasonable to expect the 
fundamental component of a resulting disturbance 
to be periodic with the approximate interval L’, 
That is, we expect that it will contain as its long 
wave-length Fourier component of appreciable 
amplitude a term with n = 2L/L*% One source of 
such characteristic maxima might be well logs that 
have been properly interpreted. Certainly this should 
be true for n because it is characteristic of vertical 
variation. It is fortunate that precise values of a 
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and B would seem to be much less important. 


USE OF STABILITY THEORY WITH LINEAR FLOW 


Symmetry makes it possible to treat many reser- 
voir problems and laboratory core flood problems 
as if the flow were one-dimensional. For this reason 
we will discuss in some detail an example long- 
core flood. The discussion shows how stability 
theory can be used to optimize solvent recovery of 
oil. It leads to results applicable to a variety of 
practical problems. Later discussion will cover the 
more complex multi-dimensional cases. 


STABILITY IN A LONG-CORE FLOOD 

The fluid and rock properties to be used in 
example calculations are given in Table 1. Density 
and viscosity depend on the composition of the 
single hydrocarbon phase. The relationships are 
conveniently given as empirical equations, Frac- 
tional solvent concentration is c, while fractional 
concentration of dissolved gas is denoted by Cy. 
By difference, the remaining material is oil. Values 
for the dispersion coefficients based on experiment 
are given in Table 1. Other materials and other 
flow conditions could require different values. 

As an example problem, we consider flow through 
a long, square-cut core about 2 in. on a side and 
positioned with a 10° dip. The core has been flooded 
with solvent for 40 days at a frontal advance rate 
of about 1 ft/day. The flood is also assumed to 
have been started by an abrupt change to solvent, 
as shown by the input curve in Fig. 1 (Curve A). If 
flow is linear and stable so that no fingering occurs, 
concentrations will be uniform across any section 
of the core. Thus, it is convenient to picture the 
result by means of a solvent concentration profile, 
as if solvent could be separately observed within 
the system. For this case, the concentration profile 
should look like Curve B in Fig. 1 after traveling 
40 ft. This profile is defined by the equation, 
x - Ugt 


2 VDit 

in which a negligible term has been dropped. ! 
The first question to be answered is whether or 

not such a flood should be stable. Using solvent 

concentration data from Fig. 1, dimensionless sta- 

bility coefficient values such as those from Eq. 3 


c = %erfc : 


. (6) 








TABLE 1 —- PROPERTIES FOR EXAMPLE CALCULATIONS 


Linear Flow 


t for © = 3.46 x 10° seconds (t = 40 days) 
D, = 0.375 uy em?/sec 

D3 = 0.0375 ug cm2/sec 

Up = 3.47 x 10°4 cm/sec (ug = 1 ft/day) 
£ = S em (k= 2 ins) 

k = 0.87 darcies 

pL = 10 exp (—3¢— 6.3¢g) cp 

Pg = (8 — 2.9€— 7.35) x 10°4 atm/em 

G =-10° 

re) = 0.25 


Two-Dimensional Flow 
As Above Except: 
t = 1.38 x 108 sec (t= 1,600 days) 
£ = 305 cm (L = 10 ft) 
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FIG. 1 — CONCENTRATION PROFILE, LONG-CORE 
EXAMPLE. 

are plotted in Fig. 2 for three different, long wave- 
length disturbances. Instability is indicated over a 
considerable concentration range in each case, as 
shown by negative values of the coefficient y. The 
results imply a substantial space region within 
which such disturbances can grow. Because natural 
core properties always are such as to initiate small 
disturbances, it is doubtful that displacement could 
have occurred as shown by Fig. 1. A more likely 
sequence would be initial instability followed by 
viscous fingering and, finally, a larger mixing zone 
than that shown. 

GRAVITY STABILIZED MINIMUM SOLVENT SLUG 

The hypothetical core flood just described can 
be modified to make use of gravity stabilization. 
Suppose conditions are maintained such that the 
function gj (Eq. 4) is positive. Then by Eq. 5, the 
concentration gradient can be made as steep as we 
like. To see what is required to accomplish this, 
Fig. 3 shows g, as a function of concentration ¢ 
for possible flood rates up. 

Since gravity stabilization requires that g, be 
positive for all c, the limiting value is obtained at 
c¢ =0 where fluid viscosity is greatest and the re- 
sponse to gravity is least. Thus, a critical rate for 
gravity stabilization is defined by setting gj =0 
and maximizing p. The result is 

unt = - 2 defd in uy” g sin 6 (7 

g 7 dc de d¢ ores ) 
For the example data, u* (Curve 4 of Fig. 3) cor- 
responds, to 0.017 ft/day. At or below this very low 
rate, complete gravity stabilization is possible 
with no reliance on dispersion. The rate is low 
enough for the dispersion mechanism to become 
molecular diffusion, but a result similar to Curve 
B of Fig. 1 would be obtained. 
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An elementary treatment of this kind of problem 
considers the entire mixing zone as a kind of ‘‘in- 
terface’’ region. This suggests as a condition for 
stability 

— A(pg) sin 6 

Au ¢ 

where A represents total change in the quantity 
across the interface. Eq. 8 is the limiting finite- 
difference form of Eq. 7 and represents a kind of 
average behavior through the interface. Used with 
the example data, Eq. 8 suggests a limiting rate 
of 0.053 ft/day. Unfortunately, correct use of Eq. 
8 requires that mixing-zone thickness be negligible. 
Because a mixed region forms almost immediately, 
only the more restrictive Eq. 7 is appropriate. 


FORMATION OF A MINIMUM-SOLVENT SLUG 

The preceding discussion showed how gravity 
could be used to stabilize displacement. Because 
such low flow rates are required, the process is 
unlikely to prove practical. The present section 
will show how to form a slug containing minimum 
solvent, yet in such a manner that flow will be 
stable at some pre-selected rate. The results apply 
to linear flow. To do this requires a somewhat 
novel approach — the use of a gradual change in 
properties of the injected fluid rather than an abrupt 
change from one fluid to another.® In addition, it is 
necessary .to consider the complete problem in which 
gas displaces solvent, which in turn recovers oil. 
Thus, the complete definitions of viscosity and 
density in Table 1 are required. The symbol ¢, is 
used for fractional gas concentration. 

The principles are as follows. Flow behavior is 
single-phase and the mixing mechanism is disper- 
sion provided (1) all fluid compositions remain in 
the miscible region, and (2) stability conditions 
are satisfied. Formation of a solvent slug which 
will satisfy these criteria involves three steps: (1) 
finding the maximum rates of change of solvent 
concentration consistent with stability for the speci- 
fied conditions and, thus, the way in which solvent- 
oil and gas-solvent transition zones are to be in- 
jected; (2) finding how much pure solvent must be 
injected between transition zones to make certain 
that dispersion mixing does not lead to immiscibility; 
and (3) finding out how much the rate can be in- 
creased as the solvent slug travels and stability 
criteria become less restrictive. These steps are 
illustrated by an example. 

Step 1 

The first step is to establish maximum concen- 
tration gradients consistent with stability. This 
can be done by applying Eq. 5 separately to solvent- 
oil and gas-solvent transition zones. The results 
specify just what the ‘‘shape’’ of the injected con- 
centration profile must be. We use the core and 
fluid properties in Table 1 and a displacement rate 
of 1 ft/day. To obtain the required data, note the 
following equivalence. 

ir ® Bc * 
oc* 1 oc (9) 
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If a concentration ct is being injected at the initial 
time t}, Eq. 10 tells what concentration ¢* must be 
injected at any subsequent time t* to keep the 
solvent slug stable but as small as possible. 


Eq. 10 is used first to construct the leading 
transition zone by which solvent is to displace oil. 
The conservative values B = 0 and n = 1 can be 
used for the perturbation properties. Eq. 4 defines 
the function g, which depends on both density and 
viscosity. To compute these functions, note that 


Cg = 0 at this leading part of the transition zone. 


It is convenient to measure time ¢* relative to a 
time ¢t¥ = 0 for which c¥ = 0. In this way we arrive 
at the solvent-oil injection profile c% shown in Fig. 
4 as a function of time. A backward time scale is 
used, listing time after start of injection. Thus, 
the profile appears the same as it would if we were 
able to observe it within the core. 

Construction of the trailing transition zone where 
gas displaces solvent is handled similarly. In this 
case Cg = ] ~C, a condition which denotes no oil 
present. Thus, the function g; has different numer- 
ical values at the trailing edge. Using Eq. 10 once 
again but with the new data, the shape of the gas- 
solvent transition zone Ct. is as shown in Fig. 5. A 
convenient starting point for the calculation is c} = 
1 at t} =0. 

The curve calculated from Eq. 10 at the trailing 
edge of the slug gives the same ¢* value for two 
different concentrations. This occurs because the 
permissible gradient becomes infinite at a concen- 
tration just below c¢=0.8. When solvent concentra- 
tions fall below this point, the mixture becomes 
gtavity stabilized. Solvent is minimized by abrupt 
termination of the slug, following the solid curve. 
The dashed extension giving double values in Fig. 
5 merely shows what the limiting behavior would be 
if the gas-solvent positions in the system were 
reversed, 
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FIG. 4 — STABLE SOLVENT-OIL TRANSITION ZONE, 
LONG-CORE FLOOD. 


Step 2 

The second step is to find a way to form a slug 
that includes the stable transition zones, and that 
will not separate into two fluid phases at some 
later stage. The method is trial-and-error solution 
of a dispersion problem. Phase behavior is assumed 
to be known as a function of c and Cg. (Note that 
the flow equations were written in terms of volume 
fractions, not mole fractions). 

We start by visualizing the two transition zones 
separated by an unknown amount of gas- and oil-free 
solvent. The transition zones are shaped as noted 
under Step 1. These zones are conveniently labelled 
as C%, for the solvent-oil transition, Cc}, for plain 
solvent and c% for the final gas-solvent transition. 
Each ¢* will be a function ¢*(7) of time 7. The time 
ty required to inject the solvent-oil transition zone 
c% is known (Fig. 4), as is the time difference ¢3 
— tg for the final gas-solvent transition cé (Fig. 5). 
The problem is to find the unknown time interval 
to—t,, during which plain solvent must be injected 
to prevent two-phase behavior at some later time. 

When these fractions of the total solvent require- 
ment are put together into a single injection pro- 
gram, the solvent concentration at any later time ¢ 
and a distance x into the system will be-given by 


aT 
tlx,t) = Pi e * (rT )Fir)dr 


t 
mm 
a J filt—to)F(r)dr 
2 





x~u,(t—t,) 
+ Kerte | = 


2,/D, (t-t,) 


x~ug(t-t) 


2,/D, \t-ty) 


ee Fe 





- %erfc 


Negligible terms have been dropped from Eq. 11. 
The functions ¢g(r) and G2(r—t2) are the transition 
zones as shown in Figs. 4 and 5. These are given 
implicitly by Eq. 10 as indicated under Step 1. The 
time variable r is a dummy integration variable, 
and the function F(r) is 
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FIG. 5 — STABLE GAS-SOLVENT TRANSITION ZONE, 
LONG-CORE FLOOD. 
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In a similar way, gas concentration is given by the 
following equation. 





c,(x,t) = % erfe 


x—Uy(t-to) 
| 


2VD, (t-t.) 


c * (tt 9)F(r)dr 


ty 
j rd 
ss tcietiaapaneae 


Oil can be obtained by difference because the sum 
of all fractional concentrations must be one. 

Starting with any reasonable estimate for the 
time tj, trial-and-error procedures will lead to a 
solution. For example, let ¢ be the time required 
for the flood to progress through the entire system. 
Using a trial value of t2, compute c and cg, for a 
number of positions x. If all ¢ and cg values cor- 
respond to a single fluid phase, the trial ¢2 over- 
estimated the amount of solvent needed. A smaller 
ta should then be used. On the other hand, the data 
could indicate that two fluid phases would be pres- 
ent at some point in the system. This shows a need 
for more solvent and a larger value of tz. The final 
time difference that proves satisfactory, tg — fj, 
will indicate the shortest time for which plain sol- 
vent can be injected between leading and trailing 
transition zones while maintaining single-phase 
behavior throughout displacement. The result, t2 - 
ty = 0, can be a correct answer for this problem and 
often may be a useful starting value for to. 

Phase behavior for the example system is indi- 
cated by the data in Table 2. These data lead to 
the injection profile of Fig. 6 (Part A) that will 
remain both stable and in a single fluid phase 
while traveling 40 ft through the example long-core 
system. This profile contains less than 8 per cent 
of a pore volume of solvent. The appearance of the 


13 





F.g. 6. Just a trace of oil in the gas-solvent transi- 
tion region is enough to reach the miscibility limit. 
Everywhere except at the single indicated point 
the mixture is safely within the single-phase region, 
even though no region of pure solvent remains. 
Step 3 

The final step is to see how much, if any, flow 
rate can be increased with time. As the flood pro- 
gresses, dispersion spreads out the solvent bank. 
Smaller concentration gradients should permit higher 
rates. The appropriate generalization of Eq. 7 for 
the critical rate becomes the minimum positive 


value of 
a k g sin 0 2p oc " op oe 
K P Oc Ox de, Ox 
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Eq. 14 includes the assumption that a = 0 and is 
valid for one-dimensional flow. It allows for the 
effect of, dispersion in stabilization, which is an 
addition to Eq. 7, and also allows for a three- 
component fluid. 

The function u* from Eq. 14 for the slug profile 
of Fig. 6 at 40 ft is plotted in Fig. 7. Because u* 
is a limit function, Eq. 14 gives the limiting rates 
(from a possibly large range of rates as shown in 
Fig. 7) for which the displacement is stable. Thus, 
the shaded areas represent rates for which the cor- 
responding part of the profile of Fig. 6 (Part B) is 
stable. As shown in Figs. 6 and 7, the critical rate 
is determined near the leading edge of the solvent 
bank. This is the smallest positive value of u* ob- 
tained from Eq. 14, as shown by the data in Fig. 7. 
Negative rates calculated from Eq. 14 have no 
meaning. For the example data, the critical rate is 
about 2 ft/day, and a gradual increase in rate through- 
out the displacement is possible. 





TABLE 2 — LIMITING COMPOSITIONS FOR MISCIBILITY, 
EXAMPLE THREE-COMPONENT SYSTEM 


profile after traveling 40 ft is shown by Part B of 








If Cg !s Se Must Be> (1-¢~¢,) Must Be < 
1.00 0 0 
0.97 0.0291 0.0009 
0.94 0.0582 0.0018 
0.92 0.0774 0.0026 
0.90 0.0966 - 0.0034 
0.89 0.1062 0.0038 
0.88 0;1157 0.0043 
0.87 0.1251 0.0049 
0.86 0. 1344 0.0056 
0.85 0.1435 0.0065 
0.84 0.1524 0.0076 
0.83 0.1611 0.0089 
0.82 0.1695 0.0105 
0.81 0.1777 0.0123 
0.80 0. 1854 ~ 0.0146 
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FIG. 6 — (A) MINIMUM SLUG INJECTION PROFILE 
AND (B) RESULTING PROFILE AT 40 FT, 
LONG-CORE FLOOD. 
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FIG. 7 — CRITICAL RATES FOR THE PROFILE OF 
FIG. 6 (PART B). 


USE OF STABILITY THEORY WITH 
TWO-DIMENSIONAL FLOW 


A mathematical solution to the problem of two- 
dimensional, stable displacement behavior has been 
obtained ! for a set of reasonable assumptions. The 
present section of this paper shows an example re- 
sult which makes use of this solution. The example 
reservoir is linear, 10-ft thick and has a 10° dip. 
The flood rate is 1 ft/day. Other properties are 
given in Table 1. While profiles are useful in de- 
picting one-dimensional flow behavior, two-dimen- 
sional flow is best illustrated as a solvent concen- 
tration ‘‘wave’’. Any cross section through the wave 
then gives a concentration profile. If fingering did 
not occur, therefore, the advancing solvent can be 
pictured as the wave shown in Fig. 8. Fig. 8 is an 
oblique view illustrating the shape of the wave after 
it travels 1,600 ft. Concentration is plotted along 
the vertical axis. The positive x-direction is to the 
right and the positive z-direction into the page; 
thus, the view is the same as if the reservoir were 
tipped on its side, For clarity, z-distances are also 
magnified relative to x-distances. Concentration 
profiles parallel to the direction of flow are given 
for distances of 0, 2.5, 5, 7.5 and 10 ft abovethe 
lower reservoir boundary. 

Similar calculations reveal that, at the short dis- 
tance of 100 ft, lateral dispersion has had almost 
no effect at all. The effect gradually increases, 
until at 6,400 ft the solvent wave is completely 
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smeared out in the (vertical) z-direction. As shown 
in Fig. 8, some segregation persists at 1,600 ft. 
The regions stable to a cos my cos mz disturbance 
are cross-hatched in Fig. 8. These exclude large 
regions of intermediate concentration with high 
concentration gradients in the x-direction. The 
steepest gradients in the z-direction occur near the 
center of the system. Here the stabilizing influence 
of gravity is felt most strongly, and the thickness 
of the unstable region is reduced to a minimum. 

It is unlikely that displacement could occur with- 
out fingering when unstable regions persist as 
shown by Fig. 8. The perturbation considered has 
a wave length of 2 ft. This is equivalent to con- 
sidering the example 10-ft section of reservoir as if 
it were made up from a stack of similar 2-ft thick 
sections. Within this stack, similar inhomogeneities 
would be (at most) 2-ft apart. Although each geologic 
formation must have its own pattern of variability, 
such values would seem to be reasonable for many 
natural oil-reservoir rocks. We conclude that, for 
stability and for maximum efficiency, the flood de- 
scribed above should have been conducted under 
more restrictive conditions. A relatively high oil 
viscosity of 10 cp was included in the assumptions 
for this example, however. A more suitable reservoir 
would contain lower viscosity oil. For example, with 
2-cp oil and no change in the other properties, dis- 
placement would be stable except during an early 
period of high concentration gradients. Provided 
initial problems could be overcome, efficient dis- 
placement would become practical. 


SIMPLIFIED METHODS FOR 
USE OF THEORETICAL RESULTS 


The principles outlined in the preceding sections 
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FIG. 8 — CONCENTRATION PROFILES AT 1,600 FT, TWO-DIMENSIONAL EXAMPLE. 





can be applied to any reservoir problem, no matter 
how complex the flow pattern may be. A more exact 
theory could have been developed by starting from 
a more detailed mathematical picture of flow within 
an imperfect system. This would add complexity, 
however, and does not at present appear justified. 
Economics will usually dictate just how detailed 
a useful theory should be. For some purposes it 
would seem that even simpler methods than those 
presented could be adopted; therefore, we will dis- 
cuss several useful simplifications. 


EQUIVALENT LINEAR SYSTEM 

Whatever the detailed pattern of a field-scale 
flood, a useful approximation is as an equivalent 
rectangular system in which flow is linear. That is, 
rows of injection wells and production wells are 
replaced by lines drawn to connect members of a 
row across the system. The distance between in- 
jection and production lines drawn in this way then 
is chosen to represent the path length along which 
dispersion mixing will occur. Thus, the mathematical 
problem is reduced in crude fashion by one dimen- 
sion, and the methods presented herein become 
applicable. Growth of the transition zone during 
stable flow is related approximately to the square 
root of the distance traveled. Thus, to minimize 
solvent requirement, a single slug should be made 
to travel the entire length of a long and narrow 
system. 


EXTENDED ‘‘RECTANGULAR’”’ SLUG 

Although the concept of the injection of a solvent 
gradient is sound, many practical problems are 
involved. These may make further consideration of 
abrupt injection of solvent worthwhile. Fingering is 
almost certain to follow such an abrupt change and 
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will increase the solvent requirement. The question 
is, how much should such an excess be? 

The first step toward a useful approximate an- 
swer is to compute the time required to inject a 
minimum-solvent transition zone, as given by Eq. 
10. In place of the gradient, however, plain solvent 
is to be injected for the same total length of time. 
Thus, the solvent injection profile becomes rectan- 
gular in shape and contains a solvent excess that 
depends in an appropriate way on the degree of in- 
stability. This makes a reasonable allowance for 
instability. Miscibility calculations for a rectangular 
slug with stable flow become particularly simple. If 
fingering ensues, however, the computed results 
could be in serious error. The excess solvent in- 
jected because of instability will be lost. A proper 
way to approximate the extended slug requirement 
is as follows. 

Step 1 

When the injected concentration profile is rectan- 
gular in shape and dispersion is the mixing mech- 
anism, solvent and gas concentration equations 
reduce to 


XU ot 
c = %erfc |—=—| 
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The first step is to compute ¢ and Cy, using Eqs. 
15 and 16 for various positions x within the transition 
zone. The time t” is the time for which plain sol- 
vent is assumed to have been injected. Trial-and- 
error procedures are used to find the smallest time 
t’ for which all portions of the solvent slug will 
remain miscible until it has reached the end of the 
system. 
Step 2 

Compute from Eq. 10 the total time required to 
inject a stable solvent-oil transition zone. Repeat 
the computation to find the time to inject a gas- 
solvent transition zone. The function gj will differ 
for the two computations. 

Step 3 

Take one-half the sum of the transition-zone times 
from Step 2 and add this to the minimum time ¢’ 
found under Step 1. The total is the estimated time 
during which solvent alone should be injected to 
form an extended rectangular slug. It includes an 
allowance for instability. 


GRAVITY SEGREGATION 


The two-dimensional result in Fig. 8 shows the 
behavior of a stable concentration wave after it has 


traveled a substantial distance. During earlier 
stages, a large amount of ‘‘extra’’ solvent appears 
necessary because of gravity segregation. For the 
example data at a distance of 6,400 ft, however, 
the profile proved nearly uniform in cross section. 
In addition, it was spread out in the direction of 
flow little more than it would have been had true 
one-dimensional flow and dispersion prevailed. 

We conclude that no allowance for gravity segre- 
gation is required when a solvent bank is to travel 
as much as 1,000 times the bed thickness. In such 
cases the problem can be reduced to one dimension. 
This does not apply to shorter, thicker systems for 
which solvent may tend to override the oil. These 
require detailed study. 


CONCLUSIONS 


A theory of stability in miscible displacement of 
oil by a solvent has been presented. Application of 
the theory leads to the definition of limiting condi- 
tions for stable miscible displacement, such as a 
minimum slug size. The criteria that must be satis- 
fied are miscibility and stability. Stability means 
that no viscous fingering will occur and that mixing 
will be caused only by dispersion. Miscibility means 
that complete recovery in the swept regions will be 
obtained. Thus, for any specified set of reservoir 
conditions an optimum use of solvent is defined. 
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The Development of Stability Theory for 
Miscible Liquid-Liquid Displacement 


RICHARD L. PERRINE 
MEMBER AIME 


ABSTRACT 


A stability theory is developed for miscible liquid- 
liquid displacement within a porous medium. In the 
usual case considered, a high-density. high-viscosity 
“oil’’ is displaced downdip by a low-density low- 
viscosity ‘‘solvent’’. Perturbation methods are used 
to find the conditions under which the spreading 
mechanism changes from the stable dispersion 
process to unstable viscous fingering. We find that 
instability is conditional, that there is a dependence 
on the shape of a disturbance leading to a‘‘diameter’’ 
effect and that very difficult experimental scaling 
problems may result. A useful consequence is the 
definition of a minimum ‘‘slug size’’ for stable mis- 
cible displacement. This should make possible op- 
timum use of the solvent process for oil recovery. 
The results may apply to many situations in which 
one fluid displaces another of somewhat different 
fluid properties within a porous medium. 


INTRODUCTION 


The process of miscible liquid-liquid displace- 
ment looked very favorable when it first receivea 
serious consideration as a means to increase pe- 
troleum recovery. Some early laboratory results 
were interpreted to mean that only a small ‘‘slug”’ 
of solvent was needed, perhaps 2 to 3 per cent of 
a hydrocarbon pore volume. This ‘‘slug’’ could re- 
cover the oil in an entire reservoir provided the 
fluids remained miscible.!»2 Only a slow growth of 
the ‘‘mixing zone’’, or gradual transition from oil 
to solvent, should occur. This would follow natur- 
ally from mixing by a dispersion mechanism such 
as that described by Scheidegger.> And, indeed, 
laboratory cores have shown this kind of behavior 
provided solvent viscosity was at least as great as 
that of the oil.24-© The same kind of behavior has 
also been observed with solvent viscosity lower 
than that of the oil after the distance traversed be- 
came large. Not all laboratory results have been 
this favorable, however.’-9 Mixing zone growth may 
become proportional to the distance traveled rather 
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than the square root of this quantity. Such behavior 
accompanies high rates in short systems of large 
diameter, provided the viscosity ratio is adverse. 
If this more rapid spreading were to persist in a 
reservoir, the solvent requirement for miscible 
“*slug’’ displacement could exceed 30 per cent of a 
pore volume. There is considerable economic im- 
portance in the difference between 3 and 30 per 
cent solvent. 

A logical conclusion is that two different kinds 
of flow behavior are possible, with each kind lead- 
ing to a different result. One kind gives efficient 
displacement, the other does not. In the efficient 
case, the solvent bank spreads out only by the mech- 
anism we have termed dispersion. Under these con- 
ditions, displacement is as near piston-like as 
possible. In the second kind of flow, found only 
with adverse viscosity ratios, viscous fingering 
occurs. That is, permeability variations cause a 
small finger of low-viscosity solvent-rich fluid to 
move ahead of its average position within the mix- 
ing zone. This creates a path of low resistance to 
flow, and an even greater amount of solvent-rich 
fluid follows. Thus, the process is autocatalytic. 
Once started, the fingering mechanism rapidly be- 
comes dominant. 

The question to be answered is this. Under what 
conditions will each of these two different kinds of 
flow occur? In particular, what conditions denote 
the transition from dispersion to viscous fingering 
as the solvent spreading mechanism? The answer 
may tell us whether or not the desirable, near piston- 
like miscible displacement process is practical. An 
answer to this problem can be obtained from theory 
by the use of perturbation methods. The procedure 
is as follows. We first formulate a mathematical 
representation of the system. Then a small disturb- 
ance (or perturbation) in solvent concentration pro- 
file is introduced and observed to see what happens. 
Unstable flow is indicated when a small disturbance 
will grow larger. This will lead to eventual viscous 
fingering. If on the other hand the disturbance dies 
out, the displacement is stable. In this latter case, 
with dispersion as the spreading mechanism, near 
piston-like displacement is possible. . 

Thus, this paper presents the development of a 
stability theory for miscible liquid-liquid displace- 
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ment. The problem differs from immiscible displace- 
ment theory !° in that no interface between fluids 
is involved. The results are important because they 
state the limiting conditions required for efficient, 
stable miscible displacement of oil by a solvent. 
Thus, they determine a minimum ‘‘slug’’ size. The 
results are optimistic to the extent that only small 
disturbances in an otherwise ideal system are con- 
sidered. A less ideal system could only lead to less 
favorable results such as the requirement of addi- 
tional solvent. When the required conditions prove 
impractical, the reservoir is not suited to the mis- 
cible ‘‘slug’’? process. An alternative recovery 
method should then be used. 


APPLICATION OF PERTURBATION METHODS 


Mathematical description of stable miscible dis- 
placement involves two transport mechanisms — 
flow, and what we will term ‘‘hydrodynamic dis- 
persion’’. 3,11-16 While such concepts as flow and 
diffusion are well understood, dispersion is rela- 
tively new. Thus, a brief description is included. 
The reader is urged to consult the references for 
thorough, rigorous development. 

We consider a system that is semi-infinite in the 
x-direction, infinite in the y-direction, and bounded 
at z°=0 and z’=L in the z-direction. The x-axis 
is rotated to a negative angle 0 below the horizontal. 
Flow is one-dimensional in the positive x-direction, 
and the respective dispersion coefficients are 
assumed constant. Fluid density is a linear func- 
tion of solvent concentration, while viscosity de- 
pendence is exponential. Although not necessary, 
these assumptions give simplicity with little loss 
in generality. The fluids are assumed incompres- 
sible, and the system homogeneous in its properties. 

We will use the following notation, where primes 
indicate actual variables and unprimed quantities 
their dimensionless equivalents. 


x =s"fL ~- ugt'/, #, 22°, 7 +9 
"fl, 4 = ugt’/L, a= u'/uo, 
= kp‘ /Lug Uo, 


= D,‘/Lug, D3 = Do'/Luy = D3'/Lug 


v'/uo, w = w' /uo, p 


Ng = ~* » ete. 
The space co-ordinates, thus, are x, y and z, and 
t is the time. The linear displacement velocity for 
stable displacement is ug, while po is the fluid 
viscosity at zero concentration and c is the frac- 
tional solvent concentration. Velocities in the x, 
y and z directions are u, v and w. These velocities 
are volumetric flow rates per unit cross-sectional 
area and per unit porosity. P is pressure, k per- 
meability, p density, g the acceleration of gravity 
and ¢ porosity. Dispersion coefficients are denoted 
by D, and letter subscripts are used to denote par- 
tial derivatives. 

To help visualize the way in which fluid moves 


through porous rock, it is convenient to think of 
the fluid as if it existed in many separate fluid 
**elements’’, each bounded by an imaginary surface. 
Thus, it is possible to discuss the motion of a 
single element much as if it moved independently 
of all others. Thinking this way and recognizing 
that porous materials are homogeneous only in a 
gtoss sense, at any time some elements will be 
moving faster than others. The relative speed will 
vary with local pore geometry. Thus, of many pos- 
sible paths leading downstream, some will prove 
to be much faster on the average than others. An 
element entering a fast path will travel further in 
a given period of time than the average of all such 
elements. But there must also be some slower than 
average paths. The net result is that a cloud of 
fluid elements starting out together will become 
spread out, or dispersed, both in the direction of 
flow and to the side. A few will end up far down- 
stream, a few will move very little and most will 
travel at nearly the average rate. This spreading 
is what is termed dispersion. 

The several theoretical treatments of dispersion 
that have been given3,11-16 lead to a common re- 
sult: to a good approximation, the dispersion pro- 
cess is mathematically described by the diffusion 
equation with constant dispersion coefficients 
(units of a diffusion coefficient). Dispersion coef- 
ficients in the direction of flow (axial) and at right 
angles (lateral) may differ substantially.15-17 Ar 
some low rate, the dispersion process necessarily 
reduces to molecular diffusion. Thus, the descrip- 
tion of stable miscible displacement requires a dif- 
fusion-type equation for solvent continuity, an 
equation for continuity of velocity and Darcy’s law 
for the equations of motion. In terms of the stationary 
dimensionless co-ordinates (x,), these become 

c¢, - Dy\c,, ~- D ~ D3c 


3° yy 22 


O wsicnnnwnnh bene 


een fy 


(kpg sin @)/ugud = 0 
winnione Ea Se 


w+ (Ho/H)p, 


navies (Ey) 


v + (Ho/u)p, = 


w + (o/p)p, + (keg cos O)/ugud = 0 


pinceeeciin. Lae 


DEVELOPMENT OF PERTURBATION EQUATIONS 
There are two kinds of displacement — stable 
and unstable — so that we are led to ask the fol- 
lowing question. Which kind of flow is more likely 
to occur, and under what conditions? Mathemati- 
cally, the problem can be stated this way. Suppose 
there is a quantity C which represents solvent con- 
centration for stable flow. It is a solution of Eqs. 
1.1 through 1.5. Suppose also that the actual sol- 
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vent concentration c differs just slightly from the 
barred quantity. Will c then change toward ¢ or 
away from it? 

This problem can be pictured for a hypothetical 
system by drawing contours along which concen- 
tration remains constant. Fig. 1 shows such a 
visualization for a two-dimensional system. Initi- 
ally, all contours parallel the z-axis in an ideal 
system, as shown in Part A. This corresponds to 
stable behavior. A more natural kind of behavior is 
shown in Part B of Fig. 1. Rather than being straight, 
the contours contain small disturbances, or pertur- 
bations, such as would result from the normal 
variation of properties within a real system. When 
instability occurs, any small disturbances (as in 
Part B) may grow and acquire a finger-like shape. 
This result is what is shown in Part C. 

Thus, the dependent variables in a real system 
(including small perturbations) are denoted by a 
set of equations of the form given by Eq. 2.18 


? 


H= ¢, By Bs O, P cnn Se 
Each barred quantity represents a stable displace- 
ment result, such as is shown in Part A of Fig. 1. 
These are solutions of Eq. 1.1 through 1.5. The 
remainder of Eq. 2 represents the perturbation. If 
Eq. 2 is also to be a solution of the equations of 
Eq. 1 with an arbitrary small ¢, after substitution 
the coefficients of «* must vanish separately for 
each n. Thus, from the coefficients of «1 = «, we 
obtain the first-order perturbation equations. 


1 1 1 1 
c, ’ ~ a hig )_D, ef = By “nr 


Ds wef 4 eu) 


+ ue, ()) + ¥e,' . 


{23 47 


+ c pw 1) == 


iimminetek 
aii 


[=e dp sin @ 
¢{[—2 = ea 
Lu» dc gp 


vionvcaeia ann a 


pie] JD = Q 


(onniin CH 


+ [ £6 dp cos @ 
jy de p 


tS cise eee 


The dispersion coefficients D, and D3 are consid- 
ered to be unaffected by the perturbation. This is 
not strictly correct because they depend on the 
displacement velocity. For cases of interest, how- 
ever, flux by dispersion is a small fraction of the 
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(C) TO DEVELOPMENT OF GROSS 
VISCOUS FINGERING 


12220? 


FIG. 1— GROWTH OF FINGERING AS SHOWN BY CON- 
TOURS. 











total. We considerably simplify the problem by this 
neglect of a second-order effect. 

So long as discussion is restricted to the first- 
order theory, the superscripts (1) for the perturbation 
variables can be dropped. The bar is retained to 
designate the properties of the stable displacement. 
At this point, we also make use of several assumed 
properties of the linear stable flow: u=1, v=w= 
0, Cy =0. Transforming to moving space co-ordinates 
x =X, —t, we obtain as the perturbation equations 

“— Dicey. st D3eyy = D3¢22 


(4,1) 


-- (4.2) 


&,c= Oncuweunataes) 


ee (io/E)P,. 


v + (Ho /E)py Dir civein ssimsisnmiennmerwh GG) 


w+ (Mo/B)ps + B36 = Onmmnnnne (4.5) 
where the functions g (of ¢) are defined by 


d inp | kg de sin & 


i ‘Cn (5.1) 


Se = 


kg decos@ (5.2) 


Se = = 
3“ uy de ¢ 
Because the set of equations of Eq. 4 are linear, 

an arbitrary initial perturbation is readily expressed 

as the sum of its Fourier components. This assumes 
that conditions of convergence and differentiability 
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are met. The equations of Eq. 4 describe the de- 
viation of a small perturbation from stable behavior. 
Higher-order perturbation theory can be applied for 
larger disturbances by also considering variables 
with superscripts of two or higher. 


SOLUTION OF THE PERTURBATION EQUATIONS 

Linearity of the equations of Eq. 4 also means 
that the solution for an arbitrary perturbation will 
be the sum of solutions for the individual Fourier 
components. Thus, we need only study separately 
the behavior of particular complex Fourier terms. 
One such follows. 


n= of] exp lif(x,y,z,t)], 


© OR Oil ccccccnnciamnainnaeee 
where 


flay. Zz, 0) = Ax + By + ANE occickheo) 

It is understood that the physical quantities repre- 
sented by Eq. 6 are real. Substituting in the equa- 
tions of Eq. 4 from Eq. 6, 


[ tf 2 DaFex ~ Dah yy - Dahan + Dif e 
+ Dy fy + D3 f2)%+e,8+2,8 =0 
AB.2) 


if + if yO + if WH Onn nnn nnnn (B62) 


Cnvcnnwnt ed) 


g, ota t iluo/Mf,p 


V+ E(Mo/ B) FyP = Oeeseneee ne nn (Bo 4) 


gs 4 B+ iluy/B) FPF Onnnnnn(B.5) 
In matrix notation, these equations are of the form, 
Ay =0, which requires either that the column vector 
¥ is zero or that |A| = 0. Thus, we obtain 


{lif, - Dy fx - D3 fyy - D3 fez + Dif? 

+Dyf2 + Df 72 +67 + FF1) 

+ Ef? +Fhf,+Gf2 +(E+G@f? =0 
nim 


where the functions E, F and G are defined as 
pe SD: 


E =_- 83 C ge teers crene 


Ss fy + Oy fy Pe 


en | | ey 


The solution /(x,y,z,t) to Eq. 9 must satisfy the 
initial condition, Eq. 7, and also the boundary con- 
ditions 


Real Cc, = 0 
} s = @, # = Leswckhl 9) 
Real w = 0 


To solve Eq. 9, we assume that the function can 
be expanded as a Taylor series. This requires that 
the complex variable / be an analytic function. If f 
is also to satisfy Eq. 7, 


G=-+ 8i C yor cn en ee ov eee © 00 00 00 00 00 
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f(x,y,z,t) = ax + by + amz + fit 
+ (QTL fS .t? + 2f2 (xx, )t 
+ 2fy2(¥-Yo)t 
+ 2foe(2-z,)tl+ x & 


vienna AEE 


The superscripts o indicate that the coefficients 
are to be evaluated at the point (9,V¥o9,2 9,0). We can 
similarly expand the known functions E, F and G of 
the stable solution. For example, 


E = E° + Eg(x-x,) + Ej(z+z,) 
+ ESt + (207! [Bo eax)? 


° 2 © 2 
+ E> 2(z-z,) + EY,t 


+ 2E°,(x—x,)(z~z,) + QE (xx, )t 


+ 2Eo ¢(Z-z,) t] a> © wot Ri 


After taking the indicated derivatives, the results 
are substituted in Eq. 9. For / to be a solution, the 
coefficients of each term in (x—x9), (x—x,)t, t?, etc., 
must vanish. This leads to an infinite sequence of 
algebraic equations. In terms of the functions, 


872 = (a2 4 B2 4 n2n2)7 1 (14.1) 
H°(x,z,t) = (a? + B2)E° + annF° 
+ (82 4 n272)G°......(14.2) 
* (a2 4 B7)E° + annF Y 
+ (B? + n2m?)G2, ete. wu0(14.3) 


The first of these equations are 


iff = - 7D, - (6? + n?n?) D, 


TE cccrininsenniion nancies 
5 Oe NE ST) 
BF g = Deninsrcnsimuenmmssminmnnsel Sy 5) 
SOL i eee |e 
iff, = - 87H} + D HE, + D3H®,) 


iS" 2{H° [20 D, + 87 2(2aE° 
+ nmF°) - 87 4(2aH°)] 
+ H3[2n7 D, + 87 2(aF° + 2n7G°) 
OE et a 
og 6 ORE oicmmmnundiedle 


SOCIETY OF PETROLEUM ENGINEERS JOURNAL 





usasatinant ss 








iawn 


Nhabes akan 













) 








8g ® Per onrninenniannimmenatliet) 


a 8 Fe iieiiwnnane A 

Before proceeding further, it is worthwhile to 
note several facts. First, if perturbation growth re- 
quires that the higher-order theory be considered, 
these higher-order terms will contribute a higher 
growth rate. For example, if a term in the first-order 
solution grows initially at the rate exp (bt), there 
will be a corresponding part of the second-order 
solution initially proportional to exp (2bt). Thus, 
any growth implies rapid growth, and it is sufficient 
to follow the behavior of the system for a short time. 
The behavior must be known at every position, how- 
ever. Because we know the functions H for every 
position (X%9, Yo, 2%), the results are conveniently 
expressed as the family of particular solutions 
{(Xo:¥o12o.t) The Taylor series then simplifies to 
a series in time ¢ only, from which we need at most 
a few terms. 


PSP get?) = ax, + By + AMZ > + fit 
+ ya" + 9" 47" 


+ * * Diigeneammmatamnnmnwannan Ree) 
The first two coefficients of t are given by Eqs. 
15.1 and 15.5. 

Thus, the first-order perturbation solution can be 
expressed in the following form, as for concentration. 


c = ¢ exp (iax + ify + innz + iy;t) 


OXP (Vat) orcccccececeereeseeecoereroeee (17.0) 


in which the ‘‘phase shift’’ coefficient y, is given 
by 
ow 
y= Real[ = ie lie oe) 


and the “‘stability’’ coefficient yp is given by 
oo 


Vp = Imaginary| = tar i~yg a"? 
R= 
In terms of the function values given in Eqs. 15.1 
and 15.5, these are 


Yr = — 4t872 {Ho [20 Dy + (2a £°+ nik? )S~? 


~ 2aH°S"4] + HS [nn Dz 

+ (aF° + 2nnG°)8~2- onnH?S74 J} 

id acai eee 
Yr = [a2 Dy + (B2 4+ n2n2)D. > S78) 

+ 4t8"?[HS + D, He, 


wun wae 


The parameters a, B and n characterize the initial 
perturbation, as do the functions E°, F° and G° the 


+ Dy H3,] 


He oe we 08 0 we we oe 00 oe oe oe oe 
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initial state of the corresponding stable displace- 
ment. 

We wish to find from Eq. 17 the conditions under 
which the perturbation becomes smaller in time. 
Because ¢ is constant and the complex exponential 
is limited in magnitude to one, stability does not 
depend on them at all. The important factor in de- 
termining growth is the exponential term that remains 
and which, in turn, depends on the stability coeffi- 
cient yp. For stability to be assured, each such 
coefficient of time in the Fourier series must be 
positive. This will guarantee that each term grows 
smaller with time. If so much as a single value of 
YR were negative, a small disturbance could grow 
at an exponential rate and eventually lead to viscous 
fingering. This means that only the least stable 
term of the original, complete Fourier series need 
be considered. If this one term with minimum yr 
gtows smaller, all other terms will do the same. 

A study of the functions involved will show that 
stability can arise from several sources, For exam- 
ple, gravitational forces can offset viscous forces 
under some conditions and provide gravity stabili- 
zation. In other cases gravity may not be quite 
sufficient, but the additional effect of dispersion 
may keep yr positive. This second case is dispersion 
stabilization. Stability is determined by the interplay 
of forces so that behavior changes with different 
displacement conditions. However, even the most 
general problem can be reduced to determining the 
sign of the coefficient of time in an exponent. 

To repeat the reasoning we have followed, sup- 
pose we are given an arbitrary initial disturbance. 
A complete description of flow behavior would re- 
quire solution of nonlinear partial differential equa- 
tions in which coefficients (such as the flow rate) 
depend on the solvent concentration answer. In fact, 
for some conditions we may even be uncertain as to 
what equations will eventually be required to de- 
scribe the behavior. But before it can die out, any 
disturbance must first become small enough for our 
linearized theory to apply. Giving this small dis- 
turbance a Fourier series representation, any un- 
stable terms with negative y would grow rather than 
die out, although there would be deviations from 
the predicted growth for moderately large amplitudes. 
The limiting conditions for stable miscible displace- 
ment can be specified after consideration of all 
possible modes of disturbance. 

It is possible to have borderline instability with- 
out viscous fingering. If growth is slow enough, a 
return to stable displacement is possible before 
any macroscopic fingers are formed. Substantial 
growth will accelerate divergence, however. Per- 
haps the true condition for prevention of viscous 
fingering in a natural system is that the disturbance 
must die out before it travels past a second pertur- 
bation source. However fingering is defined, it is 
experimentally easy to show that, once formed, 
fingers cause a radically different displacement 
mechanism.® This less efficient mechanism is like 
immiscible displacement with no interfacial ten- 






21 











sion7»8 so that, given enough solvent, almost all 
oil can be recovered. 


SOLUTION FOR STABLE DISPLACEMENT 


The perturbation parameters y and yp defined in 
the preceding section are functions of the stable 
behavior ¢ through their dependence on the functions 
E°, FS, G° and H®. Thus, we anticipate that any cri- 
teria for stability will depend strongly on the appro- 
priate kind of stable behavior. To interpret the re- 
sult implicit in Eq. 17 will require a reasonably 
general form of the solution ¢ for stable behavior. 

In most cases of stable displacement, gravity 
will be very important. Coupled with the effect of 
dispersion, gravity must prevent viscous fingering. 
The flow will then be very nearly parallel and there 
will be a nearly horizontal ‘‘mixing zone’’ between 
the solvent and the displaced fluid. This will be 
broadened from an ‘“‘interface’’ by the effects of 
dispersion. Thus, in the mathematical description 
we assume parallel flow, or flow in the x-direction 
only. Dispersion occurs in both x and z directions, 
with the behavior independent of y position. 

Any boundary condition that may be suggested 
for the inlet end of the system, x=—t, can be ques- 
tioned. We select one which in the absence of dis- 
persion would result in a horizontal interface, sta- 
tionary in a moving co-ordinate system. Thus the 
equation to be solved and its boundary conditions 
are 


Bg ORy ily By Ht Be) 


xX 22 


c(lx,z,0) = 0, x > 0; 


Mia Ftscsst) © Citic) 
go @ 
c,(x,0,t) = ¢,(#,1;¢) = O..isant20039 
For 0S tS -cot 0: 
lim clx,z>h,t) = 1; 
x~7~-t 
lim c(x,z<h,t) = 0 ..(20.4) 
; a a 
while, for t > —cot 0: 
oe Bie 8) © Recinwans (20.5) 
x—7-t- 
where 


h=1l+ttan 6 -7/2< 6<0......(20.6) 
The time ¢ is the time for which displacement at 
the constant linear velocity ug has been maintained. 

Disregarding the effects of dispersion, these 
equations describe the kind of behavior illustrated 
in Fig. 2(A) for a time t <—cot 6. A similar descrip- 
tion is given by Fig. 2(B) for times t >— cot 0. For 
the special case where 0 = — 7/2, the system re 
duces to a one-dimensional problem with step- 
function input. 

This boundary value problem can be solved by 
construction of an image system that implicitly 
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FIG. 2—STABLE DISPLACEMENT SOLUTION WITH NO 
DISPERSION—(A) ¢ <— COT 0, (B) t >—COT 


6. SHADED AREA REPRESENTS SOLVENT. 
makes use of the Green’s function for the solution!9 
ec 


and separation of variables. We construct an ‘‘ex- 
tended’’ system of mirror images of the real system. 
This implies the extended definition of gravity g 
and dip angle @ in the equations of Eq. 1 as the 
following anti-periodic functions. 


gl(z)=g, O(z)=6, 0<z < 1....(21.1) 

glz) =-g,6(z) =-0,1 <2 < 2....(21.2) 

g(z)=glz + 2n), O(z) = Olz + Qn); 
<2 Sw xk oe 


A new boundary condition for Eq. 20.1 then replaces 
Eq. 20.4 so that, for 0 <t <—cot 0, 





(E(x, h<z<2-h,t) = 1 

lim | c(x,-h<z<h,t) = 0 

sone clx,z,t) = clx,z+2n,t); 
L omet od, 2s es « 


The remaining equations apply as before. The ex- 
tended system is such that its solution automati- 
cally satisfies the conditions contained in Eq. 20.3. 
Summing the point-source contributions throughout 
the extended system and solvent ‘‘reflected’’ from 
the x =—t boundary, we obtain as the stable solu- 
tion 


= 0 
c(x,z,t) “a Z(z, t,) exp (-w2) = d¢ 





+%erfe [(x-cot 6)/2 J Di\t + cot @)] 


0 CN 
+ exp[(x+t)/D)] J Z(z, t,) exp (-d2) 3g d¢ 
+ % exp [(x + t)/D,] erfe [(x + 2t 


+ cot 0)/2,/Dy(t + cot O)] runes (22.1) 


where 
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wtp) = [x-(1-d) cot 61/2,/ Dj [t + (1-)cot 6] 
erevenmenscerecsversecesoree (22,2) 





Mo) = [x42t+(1-p)eot6]/2 VD, [t+(1-P) cot 6] 
em 
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m=— @, ee — 0, zs 9 ws + © 
‘inion nee 


When 6 = — 2/2, we obtain the important special 
case corresponding to injection of a step-function, 
for which the solution reduces to2° 


clxz,t) =% erfe ae | 
2ND, t 


+ % exp ie | erfc | =| 
Dy 24D, t 

This result can be obtained directly by use of the 
Laplace transform and arises naturally for a verti- 
cal displacement. However, it also applies when 
the injection rate is high enough for viscous forces 
to be dominant with any angle 6. In this latter case, 
dispersion must be the stabilizing mechanism, which 
may severely restrict the range of applicability. 








PROPERTIES OF THE 
PERTURBATION SOLUTION 


At this point it is convenient to note some of the 
properties of the perturbation solution and their 
relation to the real initial perturbation. Suppose 
that we had started with the real disturbance, 

c (x,y,z,0) = 2 cos ax cos fy cos n7z 
This gives eight complex exponential terms includ- 
ing all combinations of wave numbers (+ a, t 8, + n7). 
The properties of yp are such that 


Ypl+o, +6, + nm) = Yp\—, +B, -nT) 
ee ee 


Ypl(ta, +8, - nm) = ypl~a, +6, + nm) 
Similarly, we can show that for y,, 
Yrl+a +8, + nm) = ~ y,;(-a, +8, = nz) 
Yrl+a, +8, -nm)= - y;(-a, +8, + nm) 
animal ewe) 
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Thus, the complex exponentials resulting from Eq. 
23 recombine into two families of cosine terms, in 
general with different damping and phase shift. 


¢, = cos fy cos lax+nnz+y,(a,n7) t] 


exp [~yp(a,n7) t] .......00 (25.1) 


iT 


Cy = cos fy cos [ax-n7z+y (a, -n7) t] 


exp [=yp(@,-n7) t] 2. ( 25.2) 


The extended system described under ‘Solution 
for Stable Displacement”’ is also useful for study- 
ing boundary conditions on the perturbation. Take 
an initial perturbation from the cj family, cos By cos 
(ax +nmz). The part of this disturbance lying within 
the real system, 0 < z < 1, has as its mirror image 
in the right-hand adjacent strip 1 < z < 2 a member 
of the co family, given by Eq. 25.2. That is, for 
any point z; and its geometrical image z2 = 2—z}, 


cos (ax—-n72z g ) = cos [ax-n7(2-~z 1] 


= cos (ax+n7z y-2nT) 
= cos (ax+n7z ,) 


stoneninoimael aig) 


But the stable solution for ¢ (Eq. 22) has similar 
image properties as a consequence of the anti- 
periodicity of the extended system. These are re- 
flected in the functions E°, F° G° and H®°. Thus we 
can show that, when an (a,~—n7) perturbation and 
the corresponding stable solution from the 1<z<2 
strip are mapped into the 0 < z < 1 strip by reflec- 
tion, the result is identical to that for an (a, n7) per- 
turbation. This requires that the behavior of the 
superimposed systems remain the same for all time. 
We always have members of both families in any 
real solution; thus, perturbation flux into the bound- 
ary from the right is always equal in magnitude but 
oppositely directed to that from the left. This argu- 
ment is readily extended to any other ‘‘boundary”’ 
z = integer in the extended system; in this sense, 
the boundary conditions of Eq. 11 are satisfied. 

Thus, the solution to the perturbation problem 
consists of two families of perturbation waves given 
by the equations of Eq. 25. Considered in the ex- 
tended system, these satisfy the first-order pertur- 
bation equations of Eq. 3 and the initial and boundary 
conditions. The effect of a boundary is to reflect 
the wave back into the system, with an image of 
the original in the extended anti-periodic system as 
the source of the reflection. 


SUMMARY OF THEORETICAL RESULTS 


We have shown in the preceding sections that 
there are two families of mathematical solutions 
for the equations describing growth of an initially 
small disturbance. Each solution contains an ex- 
ponential factor, exp (-ypt). The coefficient of 
time yp depends on the appropriate stable displace- 
ment behavior ¢ and also on the parameters a, B 
and n of the real initial perturbation term, cos ax 
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cos By cos nnz. For a displacement to be stable, 
the coefficient yp in each such exponent must be 
positive. Thus, it is sufficient to consider the 
particular term with minimum yp. No doubt the most 
important part of this particular term is the time- 
independent element, which determines the initial 
behavior of the disturbance. This is given by 


y= D, + (A? + n?72)D, + 87 2H°) 
Cini aa 


a special case of Eq. 19.2 for which the terms are 
defined by Eqs. 6, 7 and 14. The functions used 
are to be evaluated at any point of interest (9, ¥9,Zo) 
and at the initial time with respect to perturbation 
behavior, ¢ =0. 

In terms of the stability coefficient yp, we can 
define a surface of neutral stability on which such 
a coefficient vanishes. When present, this surface 
divides the system into two parts — an exterior 
region where all yp are positive, and an interior 
region of higher concentration gradients within 
which at least one yr is negative. Thus, given an 
appropriate perturbation source, unstable fingers 
can grow from the interior region to dominate dis- 
placement behavior. For true stability, this interior 
region must be absent, in which case the surface 
YR =0 does not exist. With time, the enclosed un- 
stable concentration range will shrink. This leads 
to the first of several important generalizations. 
CONDITIONAL INSTABILITY 

The first term in yg is a positive constant rep- 
resenting the stabilizing effect of dispersion. All 
other terms depend on concentration gradients 
which must eventually become very small. This 
means that the only instability observed is a con- 
ditional instability—one which will disappear after 
a sufficiently long time. 

“WAVE LENGTH’? DEPENDENCE 

The numbers a, B and n characterize the initial 
disturbance as the sum of a series of sinusoidal 
concentration waves. One effect of the constant 
term in yr is to make long ‘‘wave lengths’’ the 
least stable modes. Thus, stability is determined 
by the ‘‘fundamental’’ Fourier mode of an arbitrary 
initial perturbation, a term with small a, B and n. 
Provided the smallest of these numbers are large 
enough, the particular perturbation can be stable 
even though displacement is not stable in the gen- 
eral sense. 

**DIAMETER”’? EFFECT 

Wave-length dependence is most important in 
laboratory-size systems. In laboratory cores, the 
longest wave length (which determines stability) is 
likely to reflect the size of the system rather than 
characteristics of the formation material. The im- 
portance of dispersion as a stabilizing mechanism 
is enhanced many-fold by such a reduction in maxi- 
mum wave length. This is shown when dimensions 
are introduced into Eq. 27 and the first terms ac- 
quire an inverse proportionality to L?, with L the 
width of the system. Experimental results showing 
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the indicated ‘‘tdiameter’’ effect have been obtained. ? 


GRAVITY STABILIZATION 

We assume a normal case in which a low-density 
low-viscosity solvent displaces oil downdip. The 
function gy (from the equations of Eq. 5) then be- 
comes positive when gravity is dominant, and nega- 
tive when gravity is subordinate to viscous forces. 
The second function of the pair g3 is always nega- 
tive. If we consider a case in which both a and n 
are positive and gy is positive, gravity effectively 
stabilizes the mode. This can be shown through the 
defining equations (Eqs. 5, 10 and 14). 

Gravity dominance will not always lead to sta- 
bility, however, as can be shown for some combin- 
ations of parameters. In particular, gravity domi- 
mance naturally leads toward instability when the 
less-dense fluid displaces updip. 


MINIMUM SLUG SIZE 

The condition for neutral stability, yp = 0, also 
defines the minimum size of a solvent ‘‘slug’’ for 
stable displacement of oil. As an example, consider 
the simplified case for which the stable flow is 
one-dimensional, Cz = 0, and a for the perturbation 
is negligible. This latter corresponds to a perturba- 
tion extending almost indefinitely in the direction 
of flow, as could follow if it neither grew nor died 
out. With all other properties specified, we wish to 
define the maximum permissible solvent concentra- 
tion gradient, cz. In magnitude this is given by 


lee) = | (8? Dg)(g,)7 4] --~..-~.(28.0) 


Maximum gradients occur at injection before any 
dispersion has occurred. Thus, Eq. 28 defines the 
steepest solvent concentration profile consistent 
with stability that can be injected under the pre- 
scribed conditions. 

Optimum use of solvent requires the smallest 
‘mixing zones’’ possible. This implies the con- 
struction of a stable solvent slug by: (1) formation 
of a transition zone from oil to solvent, such as 
would be defined by Eq. 28 for example; (2) addi- 
tion of enough pure solvent to prevent later immis- 
cibility as aresult of dispersion mixing; and, finally, 
(3) formation of a second transition zone from sol- 
vent to gas. This final transition requires re-appli- 
cation of Eq. 28. The total quantity of solvent re- 
quired can be obtained by integration over such a 
contour, thus defining a minimum slug size. 


SCALING OF LABORATORY EXPERIMENTS 
Scaling rules for laboratory miscible displace- 
ment experiments have been presented by previous 
authors.°»21 Stability phenomena impose additional 
requirements that must be met if a laboratory flood 
is to represent particular reservoir conditions, For 
stable displacement, only the dispersion and flow 
relationships between systems need be scaled. 
We expect a similar simplification for a high de- 
gree of instability. Unpublished experimental data 
tend to verify this view. Near the limits of stabil- 
ity, however, both the degree and source of insta- 
bility must be matched in scaled experiments. 
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Inclusion of the very difficult requirement of scaled 
perturbation sources is a condition which may well 
make exact scaling impossible. 
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Effect of Oil Production Rate on Performance of Wells 


Producing from More Than One Horizon 


W. TEMPELAAR-LIETZ 
MEMBER AIME 


ABSTRACT 


The performance of a two-horizon depletion-type 
reservoir produced through combination wells is 
analyzed. By introducing some simplifying approxi- 
mations, it has been possible to obtain formulas 
which are easy to handle. Only ordinary differential 
equations are used, and the development of the 
analysis can be followed without difficulty by the 
non-specialist. A rigorous analysis has been made 
of this problem. It has been found that the approxi- 
mations: introduced in the simplified analysis are 
fully justified and that errors seldom will be more 
than 2 per cent. 

The report shows the effect of the rate of with- 
drawal upon the relative depletion of the two hori- 
zons. Numerical examples are given. 


INTRODUCTION 


In its simplest form an oil reservoir consists of 
a continuous, homogeneous body of porous and per- 
meable rock, enclosed at the top and the bottom by 
impermeable material. The flow of fluid through 
such a reservoir has been the subject of numerous 
papers, and the mathematical analysis of well per- 
formance, pressure build-up, etc., has reached a 
high degree of refinement. One of the basic assump- 
tions introduced into these analyses is the homo- 
geneity of the reservoir. In non-homogeneous reser- 
voirs, these computations still can be carried out 
satisfactorily, provided small irregularities in the 
physical properties of the rock are distributed in a 
random fashion. 

The flow equations are derived as solutions of 





Paper originally published in 1953 as Shell Oil Co. report. At 
that time a search of the literature revealed no mathematical 
analysis of the production from a multiple-zone reservoir as 
herein described, Since then, the report has been referenced in 
later studies on related subjects, 

At the 34th Annual Fall Meeting of SPE in Dallas in 1960, 
the paper by C. S. Matthews, etal (page 43 of this issue of 
Soc. Pet. Eng. Jour.) was presented. A portion of that paper is 
devoted to a rather detailed analysis of the report, and it was 
felt that publication of the 1953 report might be desirable. 

Consequently, the author submitted this paper—essentially in 
its original 1953 form — for presentation at the Regional SPE 
Meeting in Pasadena, Calif. on Oct. 22-23, 1959. Subsequent 
constructive comments by the Society’s Transactions Editorial 
Committee are gratefully acknowledged and have led to modifi- 
cation of the text to its present form, 
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the continuity equation, and their solutions are 
dependent upon the assumption of many constant 
factors in the reservoir: permeability to oil, com- 
pressibility of the reservoir fluid, viscosity satura- 
tion, etc. When these are assumed to be constant, 
the continuity equation can be solved rigorously. 
When these assumptions do not hold and corrections 
are introduced, it becomes very difficult and often 
impossible to solve analytically. 

The present analysis is concerned with production 
from a two-layer depletion-type reservoir in which 
the layers have different permeabilities. Instead of 
using the continuity equation, results were derived 
from some simple physical ideas. In a single-layer 
reservoir in which the oil is at a pressure above 
the bubble point, the rate of production often has 
been assumed to be directly proportional to the 
pressure drawdown. In the simplified analysis it 
has also been assumed that the drop in average 
reservoir pressure is directly proportional with the 
cumulative withdrawal. These considerations in 
the past have been the basis of the analysis of 
reservoir performance and have been applied to a 
variety of problems, often with marked success. In 
the problem under investigation these relationships 
are shown to be compatible with the equation of 
continuity, but the approach has the advantage of 
leading to ordinary differential equations in which 
changes in constants can be introduced without 
making the analytical solution impossible. 

It has the additional advantage that it can be 
understood more readily by the large group of engi- 
neers who are not reservoir specialists. For that 
reason, detailed steps in the mathematical devel- 
opments have been shown in full. 


SCOPE OF ANALYSIS 


The analysis in this paper has been applied to a 
two-layer depletion-type reservoir in which both 
horizons are produced simultaneously in the same 
well. Because of the difference in permeability, it 
is obvious that the depletion of the two zones will 
proceed at different rates. For various problems it 
may be of interest to know the degree of depletion 
of the two zones at any arbitrary time and to in- 
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vestigate the possibility of affecting this relative 
depletion by operating practices. One of the main 
factors under control of the operator is the rate of 
withdrawal, and the effect of various rates upon 
the distribution of depletions is demonstrated. In 
the development of the analysis, some basic assump- 
tions have been introduced which now will be 
discussed. 


BASIC ASSUMPTIONS 


The following assumptions have been made. 

1. Each of the two layers producing simultaneously 
in one well has an equal (unit) thickness. The rea- 
soning followed in the computations should not 
change when layers of different thicknesses are 
assumed, 

2. The two layers have the same initial reservoir 
pressure. 

3. The P-V-T data of the fluids in the two reser- 
voirs are equal. 

4, The well under consideration is one in a regu- 
lar drainage pattern, and surrounding wells are pro- 
ducing at comparable rates. 

5. The structure is flat, and gravity effects can 
be ignored so that radial flow formulas can be applied. 

6. The reservoir is the depletion type. 

7. The present analysis is mainly concerned with 
the period in which the producing GOR is equal or 
close to the original solution GOR, up to the time 
that free gas starts moving; and, experience has 
shown that during that period the relationship be- 
tween pressure and cumulative production can be 
approximated closely by a straight line. In this 
analysis, it has been assumed that such a straight- 
line relationship does exist. 

8, Although the permeabilities of the two layers 
may differ widely, it often is found that the porosities 
show a random scattering. It has been assumed for 
this analysis that the average porosities of the 
layers are equal. 

This assumption will lead to a simplification of 
the formulas. For the case of different porosities, 
adjustments in the analysis easily can be made. 

9. The instantaneous rate of production is directly 
proportionate with the difference between the aver- 
age pressure of the reservoir and the producing 
pressure in the well. This assumption will be dis- 
cussed in detail in the Appendix. 

10. The drainage reservoirs of the individual 
wells are approximated by circles with radii equal 
to one-half the distance between the wells. 


MATHEMATICAL ANALYSIS 


Let it be assumed that a well is producing from 
two zones simultaneously, and let the two zones be 
referred to as Zone I and Zone II. 

The recovery per psi pressure drop for each of 
the zones will be represented by A barrels (Assump- 
tions 7 and 8). 

The proportionality factor according to Assumption 
9 will be referred to as the productivity index and 
will be represented by i, for Zone I and iz for Zone 
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Il. The nature of these i-values will be discussed 


further in the Appendix. 

The instantaneous rate of production from the 
two zones combined will be denoted by q B/D, and 
in this analysis it has been assumed that this rate 
is constant throughout the period under considera- 
tion. 

The original reservoir pressure is represented by 
P,; and at a time ¢ the average reservoir pressure 
in Zone I will be p,; the average pressure in Zone 
II will be represented by p2, and the back pressure 
in the producing well by p,,, all expressed in pounds 
per square inch. 

At that time the cumulative production from Zone 


I will be equal to 0, = A(P, - P ,)- 


Similarly, for Zone II, = 
Q, = AE, — Po). 
Therefore, the cumulative recovery from both zones 
at time ¢ can be represented by 
O=Q,+Q,=2AP, - A (p, + pp). 

Since it has been assumed that the well has produced 
at a constant rate of g B/D, 

2AP, — Alp, + Pod=qt oe eee eee ee (1) 
The rate of production from Zone I at the time ¢ will 
be equal to 








Wy! 
—! = i, (B, - b,,), 
dt 
and because of Assumptions 6 and 8 
Qa - dp 
— = i, (Py — By = Ae ee ee eee (2) 
dt dt 
Similarly, for Zone II, 
dQ = dp, 
2. : a 
—= = 1, (Py — by) = -A — ha hag araman aan (3) 
dt dt 
Therefore, 
dQ, dQ 
Bind Pie ee 
dt dt at 
or “ - 
AP, + Pp — (+ RY = 4: 
so that 
gee: xe Aneedoaee (4) 
Ww 
4 + ty 
Substituting Eq. 4 in Eq. 2, 
ae i aS le 
At = isp, ~ + — (ib, + 4P2- +--+ - (5) 
dt i +h 
From Eq. 1, 
2P,-(p,+pJ=4t, 
or A 
= — 4 ; 
b= 2% -h oe bela eA elete maa rees (6) 
A 


Substituting Eq. 6 in Eq. 5, 























ii oe .t 1 
ce GR #6, = Lh 4t@, 
i, +i, A ith 
or 
dp 24, 4 1, 
14 12 p,+ 142 t+ 
dt A(i,+i,) Ai +i) 
qi,-21,i,P 
Peis a a ee (7) 
A(i, +4) 
Let 
24% 
a * =C, 
A (i, + t) 
pg 
75% me 6 
A (i, + i,) 
“eee 
dies | *s > ‘o =C, 
A (i, + i) 
Then, 
dp, Cok ee 
— 1+ gi + C,=0 ae ee ee a ee (8) 
t 


This equation is linear and can be solved as fol- 
lows. Assume py = uv. 
Then, dp dv du 
Eh Te | eS | EO Re Oe ee (9) 
dt dt dt 


Substituting Eq. 9 in Eq. 8, 


dv du 
u—_+v——+C.uv+C, t+C, =0, 
1 2 3 
- at at 
dv du 
u(—+ Cv) + (v__+ C,t + C,) 2505.4. 10) 
at dt 


One now can assume one relationship between u 
and v. dv 


M—t Cw) © 0. 6 eer edn (11) 
dt 
Eq. 10 now can be written 
du 
y—+Ct+C, = 0....... (12) 
From Eq. 11, . 
a “Cv 
dt 
-C,t 
ak emma ren See eer ta (13) 


Substitution of Eq. 13 into Eq. 12 gives 


-C,t du 
e€ —+Ct+C, = 0, 


dt 


so that 


Gai 
2 * «¢ 


c C.t 
gam? (C,t-le } 


das eatt (14) 
2 
Cy C, 


4’ 


in which C4 = constant of integration. 
With py =uv and v=e7Ct! Eq. 14 can be trans- 
formed to 


oe a C.t 
pb, =-—(C,t-1) - + + Cye sas Pe 


2 
Ci C, 


Fort=o is p; =P, or 


Cc G 
Ci = P - aah + wat - 
4 (e) 2 C 
Cy 1 
Substituting the values for Cy, Cg and C3, 
q (i, —#,) 
CLO mle (hI eae KER (16) 
4 ts i, . 
Substituting the value of p,; in Eq. 6, 
a = q 
P, = 2 a - p, - —t 
A 


using the proper values, 
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Eq. 17 gives the value of the average pressure 
in Zone II = po at any given time ¢, and substitution 
of this value into Eq. 6 gives the corresponding 
value for pj. 


p, =P ." re 
A isa eee ee, 
2A 4i, ty 


NUMERICAL EXAMPLES 


To illustrate the theory, numerical examples 
have been worked out. The following basic data 
were used. Porosity of reservoir rock = 0.200; water 
saturation S,, = 0.300; oil saturation S$, = 0.700; 
compressibility of water C,, = 2.5 x 10°® vol/vol/psi; 
and compressibility of oil Cy = 10 x 10-® vol/vol/ 
pSi. 

Therefore, effective compressibility referred to 
oF Go * Pa x Sw = 11.09 x 10-© vol/vol/psi. 


Oo 
In the first example it is assumed that the original 


reservoir pressure was P, = 6,000 psi while the 
bubble point of the oil was Py, = 3,000 psi. Further 
assumptions include: zonal thickness ) = 100 ft; 
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drainage radius rp, = 500 ft; well radius r,, = 0.25 
ft; permeability Zone I ky = 100 md; and permea- 
bility Zone II ka = 10 md. 

Considering the period in which the reservoir 
pressure is still above 3,000 psi, Fig. 1 was con- 
structed, For this plot, the value of 11.09 x 10-6 
vol/vol/psi was used, It will be observed that, while 
in the early stages the rate of production from Zone 
I represents the bulk of the well’s production, grad- 
ually the production from Zone II increases while 
that from Zone I decreases. 

After having produced about 30,000 bbl from the 
well, each zone is contributing half the production 
of the well. It also will be observed that from that 
time the lines run practically parallel. Theoreti- 
cally, the lines should come together after an in- 
finite time. At the time Zone I reaches the bubble 
point, it will have produced a cumulative production 
of 65,000 bbl. If the well had been producing at a 
rate of 100 B/D, Zone II would have produced about 
63,000 bbl, while at a rate of 500 B/D this recovery 
would be only 55,000 bbl. 

To demonstrate the effect of changes in one of 
the controlling factors, an expansion factor of 50 x 
10-6 vol/vol/psi has been used. From about 100 
P-V-T analyses carried out by Shell, this figure 
represents a good average value for the expansion 
of oil in the first period after it passes the bubble 
point. 

Fig. 2, therefore, represents the performance of 
an oil reservoir that has the same properties as the 
one considered in the first case, with the excep- 
tion that the oil has a bubble point equal to the 
original reservoir pressure. Again considering the 
time that 65,000 bbl have been withdrawn from 
Zone I (maximum expansion oil of Case 1) it now 
will be seen that, if the well had been producing at 
a rate of 500 B/D, Zone II would have produced 
27,000 bbl or less than half of Zone I (Point Z of 
Fig. 2). 

No corrections have been made for the gradual 
change in relative permeability because both are 
changing downward and the magnitude of these 
changes for 100 and 10 md may affect the difference 
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in depletion one way or the other. In those cases 
where relative permeability data are available, 
corrections in the permeability can be made if con- 
sidered necessary. 


GENERAL REMARKS 


The analysis carried out in this paper essentially 
is based on steady-state flow. As such, it should 
find application in depletion-type reservoirs, i.e., 
in reservoirs with closed boundaries and in which 
no edge water encroaches. To verify the validity 
of this approach for the problem under consideration 
in which the relative rates of production from the 
two zones are controlled by continuously changing 
pressure conditions, the same problem has been 
solved rigorously. It has been found that, with the 
exception of the first few hours of production, errors 
of less than 2 per cent are introduced by the steady- 
state approach. 


CONCLUSIONS 


1, An analysis has been made of the performance 
of a depletion-type reservoir consisting of two hor- 
izons of different permeabilities, producing through 
wells penetrating both horizons. 

2. It has been found that simple formulas, based 
on steady-state conditions, lead to answers that 
are within a few per cent from the results of com- 
putations based upon transient conditions. 

3. It is demonstrated that the relative depletion 
of the two horizons can be affected to a large ex- 
tent by adjusting the rate of production from the 
well, 


NOMENCLATURE 


A_ = recovery in bbl per psi pressure drop for each 
zone 

Cw = compressibility of water, vol/vol/psi 

Co = compressibility of oil, vol/vol/psi 

































































f = porosity of reservoir rock, fraction of bulk 
volume 
h == zonal thickness, ft 
iy = productivity index Zone I 
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FIG. 2—CUMULATIVE PRODUCTION FROM ZONE II 
VS CUMULATIVE PRODUCTION FROM ZONEI. 
A = 108.5 BBL/PSI; #1 =5.17 B/D-PSI; i2=.517 
B/D-PSI. 





12 = productivity index Zone II 
k = permeability, md 
ky = permeability, Zone I, md 
ko = permeability, Zone II, md 
P, = original reservoir pressure, psi 
Py = bubble point of oil, psi 
py = average reservoir pressure Zone I, psi 
p2 = average reservoir pressure Zone II, psi 
Pw = pressure in production well, psi 
Q = cumulative production from well, bbl 
Q, = cumulative production Zone I, bbl 
cumulative production Zone II, bbl 
= rate of production of well, B/D 
= drainage radius, ft 
= radius borehole, ft 
= water saturation, fraction of porosity 
= oil saturation, fraction of porosity 
= time 
viscosity of fluid 
density of fluid 
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APPENDIX 


In the preceding analysis, some assumptions 
have been made which call for further investigation. 
In the development of the formulas, it has been 
assumed that the productivity-indexes 7, and iz can 
be considered as constant multiplication factors. 

It will be shown first that, if the flow pattern at 
any time can be considered as one of steady state, 
this assumption is justified. 

Under steady-state conditions the pressures 
throughout the drainage-area are falling equally; i.e., 
“be is not a function of position. 

The density of the fluid can be expressed as 


p=pe ~c(P —p) 


for p < P,. The initial pressure and density will be 
taken for the reference state (Pj, po). The com- 
pressibility c is assumed constant and sufficiently 
small for the approximation: 


p=p,ll1-c(P ~p)l. 
In the steady state, the mass rate at the wellbore 
face can be expressed by close approximation as 


ant op (2 . 


in which rp = boundary radius, in this case the 
drainage radius of the well, / = porosity of the rock 
and b = thickness of the zone. From Eq. 19 one 
finds dp dp 

—= p,c— 

dt dt 
Therefore, . dp 
m= nr bf (cp, —), 

dt 


or Z 

— q= mrthfc, 

b 
. t 

in which q is now expressed as a volume rate and 
equal to the rate in volumes of original oil per unit 
time. 

The rate past some circular boundary of arbitrary 
radius 7 will be — 
(4 -r* ) 
The rate of production past a radius r can be ex- 
pressed, according to Darcy’s law, as 

) 
es i i Pe 22). 
12 ph Or 

in which k is permeability and p is viscosity, so 
that OP qe 1 ” r ' 
Or 2a kb r “ 
By integration the following expression is obtained 
for P,. 


forr=r,, P, = 


qu 

P = (In 7 

Omkh Oe 
Subtracting Eq. 25 from Eq. 24, 


2 

ne ea ae 
27kh om 

Introducing a volume average pressure P, for any 

radius r defined by 1 


V, 


r £ 
w 


in which V, is the volume of the drainage area from 
the wellbore to the radius 7, the following equation 
can be written. 


Pe 





2 


m(r?— 12) bf )2m 2r} 


1 Tp : r pte r? 
f (In _ ~___ ) (2 rhfdr) 
kh 2 r 


r 


+ J ©. (2nrhfdr) 


w 


r2 
2nkbh ) 12-1? 
. w * 
Substituting r = rp and realizing that under normal 


conditions od K rf, this formula can be written with 
close approximation as 


- igs x 3), 


P=p ao 
‘ i 2nkh ‘ 4 


Comparing 
7 2akh 


SOCIETY OF PETROLEUM ENGINEERS JOURNAL 





with Eq. 26, 


_ Gu 
P« i= In 


b 2nkh 


it appears that the difference between using P,, and 


— %*. . 
P,, when substituting normal values for -b is of 
secondary importance. 


Eq. 31 may now be written 
q-i(P, — by)» 


in which 


2nkh 


It appears that the value of i is not a function of 
time and, therefore, is constant for the reservoir so 
long as k and p do not change. 

The mass in the reservoir at any time ¢ can be 
derived from Eq. 19. 


r r, 


b b if 
J pav=o, f dV-C p, pf dV + 
r r 


T 
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Using the definition of Eq. 27 for P,, Eq. 35 can be 

rewritten. Let m represent the mass in the reservoir 

at time ¢ and m, the original mass; then, 
m=m,—Cp,F, | * © Pot 


—-—m —_ 
=CV, (P,-F), 


Po 
and with 
V,, = Mb OG ~ 40 )- O = Cabf (—7,) (Fy ~ Fs 


Q- A(P,~F ) 


in which Q = cumulative production in volumes of 
original reservoir oil and A is a constant. P,, is 
the same pressure appearing in Eq. 33. 

Hence, if we can assume that, if the formulas 
developed for steady-state conditions can be applied 
to the problem under consideration, the following 


main equations can be written. 
Q= A(P, -P) 


dQ 
dt 


These equations are essentially the same as those 
used in the main body of the report. 


i = - 


Ww 


q 





Linear Water Flood with Gravity and Capillary Effects 


S. A. HOVANESSIAN 
MEMBER AIME 

F. J. FAYERS* 
JUNIOR MEMBER AIME 


ABSTRACT 


The one-dimensional displacement equation for a 
homogeneous porous medium, including the effects 
of gravity and capillary forces, has been solved by 
a numerical method, A finite-difference scheme is 
developed for obtaining saturation, pressure and 
fractional flow profiles in waterflood recovery prob- 
lems. From the numerical examples given, it is 
concluded that the gravitational forces have a pro- 
nounced effect on the saturation profiles and the 
pressure distribution curves of the system. 


INTRODUCTION 


Within the past 20 years, a number of papers have 
appeared in the literature dealing with the quantita- 
tive treatment of waterflood recovery problems. In 
their celebrated paper, Buckley and Leverett! de- 
scribed a method for calculating saturation profiles 
when the effects of capillary pressure and gravity 
are excluded. Terwilliger, et al,? included the effect 
of gravity in their theoretical and experimental in- 
vestigation of oil recovery problems and obtained 
close correlation between experiment and theory. 
Welge? described a simplified method for obtaining 
the average saturation and the oil recovery from an 
oil reservoir. The effect of gravity can be included 
in this calculation. 

More recently (1958) Douglas, et al,4 presented a 
method for calculating saturation profiles which in- 
cludes the effect of capillary pressure. The authors 
start with the one-dimensional displacement equation 
(which is nonlinear in the derivative) and, by a 
change of variable, transform this equation to a 
semi-linear partial differential equation. This equa- 
tion is then solved by a finite-difference method on 
a high-speed digital computer. Fayers and Sheldon5 
obtained the solution of the one-dimensional dis- 
placement equation by difectly replacing the differ- 
ential equation by a finite-difference equation. The 
effects of capillary pressure and gravity were in- 
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cluded in their solution. Here, since the movement 
of the foot of the saturation front is governed by a 
separate equation, the elapsed time in attaining a 
particular saturation profile cannot be obtained. 
The results of these two papers indicate that the 
inclusion of capillary forces eliminates the triple- 
valued Buckley-Leverett saturation profiles. 

In the present paper, the one-dimensional dis- 
placement equation for a homogeneous permeable 
medium, including the effects of capillary and gravity 
forces, is solved. 

The method of solution resembles that described 
in the paper by Douglas, et al, since it involves a 
similar transformation of variable. However, it ex- 
tends their work in that the effect of gravity is in- 
cluded and, in addition, pressure profiles may be 
calculated. Also, the functions of saturation [k,,(S), 
kyyfS) and P.(S)] required are entered in tabular 
form rather than as low-order polynomials. This 
simplifies data preparation and permits the use of 
a greater range of function types. 

Unlike the method of Fayers and Sheldon, the 
corresponding elapsed time required for the devel- 
opment of each saturation profile is calculated and, 
also, saturation profiles may be calculated after 
breakthrough. 


THEORETICAL CONSIDERATIONS 


The dimensionless form of the one-dimensional 
displacement equation for a homogeneous porous 
medium as obtained from Darcy’s law and material 
balance relations can be written as follows.> 


as , dG(S) OS C) | 3 


+N: —— 16S) —T=0 
aT dS ax , 


° ax dx 

ep 
where S = water saturation, ¢ 
T = dimensionless time = = . 


G(S) = gravity function = F,,(1 — Ngkyo), 


F ,, = fractional water flow = ( 


kew Bo 


Ng = gravity constant = [(py— po) gk sin 0)/ 


Hod 
X = dimensionless distance = x/L, 
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No. = capillary pressure constant = 





C(S) = capillary function = F,,k 


Pe = Pw — Po, and 
fw = actual fractional water flow = G(S) - 


Os 
NoC(S)-2 


The solution of Eq. 1 will give the values of satur- 
ation S as a function of time T and position X. This 
equation is a second-order nonlinear partial differ- 
ential equation and it cannot be solved by the pres- 
ent classical techniques. However, it is possible 
to remove the ‘‘nonlinearity in the derivative’’ from 
the equation by introducing the following transfor- 
mation. 4,6 Ss 


riS)= > C(S).dS. 


Si 

ye at a » (2) 

where tuff ** c(S)dS, 
S 


i 
S; = initial water saturation, and 


Sor = residual water saturation. 
The introduction of transformation Eq. 2 into Eq. 


1 yields . 
1 Or . 1! AG(r) dr OT 
inte Gas Geis oe ae 8 
Cir) OT Z dr ox i ax2 9 
se 


In Eq. 3, the capillary pressure and gravity func- 
tions are both expressed as functions of the trans- 
formed variable r. This equation is nonlinear but 
not in the derivative. 

The partial differential equation (Eq. 3) can be 
approximated by the finite-difference equation, 7*8 


, s Ss 
a Tr -Tn = A(r?) [ate 





AT 24x 





- | 
— (Ax)? 
TST ete Te Se eee ee s+ 
The space subscript is n. Primed and unprimed 
quantities denote, respectively, values at the begin- 
ning and end of the time step AT. Also, 


eri + ’) aAT 
=F, 
ll oT 





« 
Tr = Om +(I1-a) tp 
~« 6) 


where a = time centering parameter, 0 <a<¢ 1 (usu- 
ally set at 0.5) and 


C(r*) WG(r*) 
z dr* 





A(r*) = and B(r*) = No C(r*) 


Using Eq. 4 to approximate (07, /0T)’ in Eq. 5, Eqs. 


4 and 5 combine to give 
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(Ax)? 2Ax (Ax)2 
as A(rx 7 
+ adT + Trot = ™ 
(Ax)? 2Ax 


ee Se a ee a - (7) 
Denoting the coefficients of 7,-1,7, and 7,44 by 
Gn, by and cy, respectively, and setting d, = r4, 
we can write Eq. 7 as 
An Tr-1 t+Ontn +Cn tng) = dp, Nel, n 

The boundary conditions are fixed as follows: 
(1) the fractional flow of water /, is always unity 
at the inflow boundary; and (2) the saturation at 
the outflow boundary must build up to a value Sout 
before water will flow out of the system, thereafter 
being calculated. 

The solution of the finite-difference equations 
(Eq. 8), together with Eq. 6, will give the values 
of r as a function of time T and distance X. These 
values are in turn transformed back to saturation 
values by means of Eq. 2. Since the numerical pro- 
cedure involves a great number of algebraic calcu- 
lations, an electronic digital computer was used 
for obtaining solutions. Equations also are written 
for obtaining the fractional water flow-.and the pres- 
sure distribution across the system. 


DISCUSSION OF RESULTS 


Saturation, fractional-flow and pressure-distribu- 
tion curves are obtained for oil-wet and for water- 
wet systems, with angle of inclination and injection 
flow rate as parameters. The relative permeability 
and capillary pressure data of the two systems are 
the same as given in Ref. 5. 

Figs. 1 through 3 illustrate the effect of gravity 
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on the saturation, pressure and fractional flow pro- 
files of an oil-wet system. From Fig. 1 it is noted 
that, for a constant injection rate g = 10-4, the sat- 
uration profiles move slower for an upward angle of 
inclination (6 = 60°) than for a horizontal or declined 
system. The areas under the saturation curves for 
the same value of T are equal, since T represents 
the pore volume of water injected. The initial satur- 
ation line starts from maximum saturation (1-S,,) at 
the input face and drops to a value slightly greater 
than the connate-water saturation and maintains 
this value for the remainder of the system.* Fig. 2 
gives pressure distribution curves for the same 
system. It may be seen from this figure that the 
greatest pressure drop occurs for an inclination 
angle of +60°. Note the slight dip in the pressure 
curves at a distance of about 0.4 and 0.8 for T = 
0.2 and 0.4, respectively, corresponding to the po- 
sition of the saturation fronts at the times indicated 
(Fig. 1). Note also that almost no pressure drop 
occurs for a dipping angle of 60° below the hori- 
zontal, indicating that the pressure difference re- 
quired to maintain the specified flow rate is approx- 
imately offset by gravity forces. Fig..3 gives the 
actual fractional flow-vs-saturation curves for the 
cases illustrated by the previous figures. From Fig. 
3, using Welge’s method of calculation, we observe 
that the breakthrough saturation increases as the 
dip angle increases (from —60° below horizontal to 
+60° above). This is intuitively reasonable because, 
in general, the greater the slope of the test section, 


the greater the tendency for the injected water to 
accumulate and, thus, to increase the water satura- 
tion rather than to move the saturation front forward. 





*The finite-difference scheme employed requires that all 
saturation values be at least slightly greater than the connate- 
water saturation. This can be seen from Eq. 4, which represents 
the rate of change of r(r = saturation), This rate depends on the 
values of A(r*,) and B(r => which, in turn, depend on the value 
of r at the beginning of the time step and, hence, are zero for 
connate-water saturation. 
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Figs. 4 and 5 give the saturation and pressure 
profiles for a horizontal section at two flow rates. 
Comparison of Figs. 4 and 1 reveals that, as the 
flow rate increases, the breakthrough saturation 
increases. This effect has also been observed by 
Perkins9 in his experimental investigations. Fig. 5 
represents the pressure distribution across the test 
section for the same input conditions illustrated in 
Fig. 4. Observe that the greatest pressure drop 
occurs for the highest flow rate, and note the dips 
on the two pressure curves corresponding to the 
flow rate of g = 1072. It may be observed (Fig. 4) 
that these dips occur at the location of the satura- 
tion front at the times indicated. 

Figs. 6 and 7 represent the saturation and pres- 
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WET SYSTEM. 
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sure profiles for a water-wet system. From Fig. 6 it 
is noted that the saturation profile for the higher 
flow rate (dashed liné) is steeper than that of the 
lower flow rate. Comparing Figs. 6 and 4, note that 
the saturation distribution lines of the former have 
a greater slope than those of the latter. This comes 
mainly from the capillary pressure curve because 
the capillary pressure curves of water-wet and oil- 
wet systems differ greatly, while the relative per- 
meability curves of the two systems have substan- 
tially the same shape. Fig. 7 represents the pres- 
sure distribution across the system. Here again, 
the higher pressure drop occurs in the case of higher 
flow rate, and the dips in the pressure curves cor- 
respond to the location of saturation fronts. The 
effects of gravity forces were investigated for 
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water-wet systems and no significant change was 
observed in the saturation, pressure or fractional- 
flow profiles. 

In obtaining these solutions, AX was set at 0.04 
and AT was set at 0.01. The material balance error, 
using these space and time increments, was less 
than 1 per cent. 


CONCLUSIONS 


From the results of this study, the following con- 
clusions can be made. 

1. The inclusion of capillarity eliminates the 
triple-valued Buckley-Leverett saturation profiles. 

2. The inclusion of gravity significantly affects 
the rate of movement of saturation profiles for a 
certain range of injection flow rates. 

3. Gravity forces have a pronounced effect on 
pressure distribution curves for a certain range of 
injection flow rates. 

4. Injection flow rates have a pronounced effect 
on saturation profiles and pressure distribution 
curves. 

5. The initial saturation distribution had little 
effect on the computed results [S(X), P(X), f,(X)] 
when the initial water in place was small compared 
to the water in place at the time under consideration. 


NOMENCLATURE 


centering parameter 
= capillary function 
fractional water flow 
actual fractional water flow 
= gravity function 
g = gfavitational constant (cm/sec?) 
= permeability to oil and water (darcy) 
= relative permeability to oil and water 
absolute permeability (darcy) 
length of system (cm) 
= capillary pressure constant 
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= gravity constant 

= pressure in oil and water phases (dyne/ 
cm2) 

= capillary pressure normalization constant 
(dyne/cm?) 

= dimensionless capillary pressure 

= volumetric oil and water flow rates (cc/ 
cm2)/(sec) 

= transformed variable 

= water saturation 

initial water saturation 

= residual oil saturation 

= saturation at the outflow boundary 

= time (sec) 

= dimensionless time (pore volume of water 
injected) 

= distance (cm) 

= dimensionless distance 

= constant of transformation 


1-So, 
f C(S)ds 
Sj 
= viscosity of oil and water (cp) 
angle of inclination 
= porosity 
= density of oil and water (gm/cc) 
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Thermal Conductivities of Porous Rocks Filled 
with Stagnant Fluid 


D. KUNI! 


J. M. SMITH 


ABSTRACT 


Effective thermal conductivities of sandstones 
filled with stagnant fluids were measured using a 
steady-state technique. Data were obtained for 
seven sandstone samples, taken from four different 
locations and ranging in permeability from 18 to 
590 md. The measurements with gases (helium, 
nitrogen, air and carbon dioxide) covered a pressure 
range from 0.039 psia to 400 psig. Data were taken 
for four liquids — n-heptane, methyl alcohol, 79.8 
weight per cent glycerol-water solution and pure 
water at atmospheric pressure. 

The experimental results were used to evaluate 
the theoretical equations for predicting stagnant 
conductivities developed earlier.4 The low-pressure 
measurements permitted evaluation of the consoli- 
dation parameter hpDpy/kg (necessary to utilize the 
theory) for the various types of sandstones. Using 
these characteristic values, the theoretical equa- 
tions correlated well with the experimental conduc- 
tivity data for the several fluids and rock samples. 


INTRODUCTION 


An aspect of heat transfer in solid-fluid systems 
of considerable current interest is the effective 
thermal conductivity of porous media, The stimulus 
for study of the subject arises from the need for 
sound procedures for designing thermal methods of 
petroleum production. The general system occurs 
when there exists a flow of fluid through the pores 
of the solid material. However, a logical starting 
point in developing a theory for predicting the ef- 
fective thermal conductivity in the general system 
is to attack the special case when the porous solid 
is filled with stagnant fluid. Since the flow rates 
anticipated in thermal production processes are 
very low, such stagnant conductivities k@ are also 
of practical significance. 
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Kunii and Smith4 recently reviewed the existing 
data for stagnant conductivities and proposed a 
theory for heat transfer in porous rocks. This leads 
to equations for predicting the stagnant conductiv- 
ity as a function of the properties of the fluid and 
solid phases. The existing experimental informa- 
tion!-3,5,7,8 covered a range of porosities and 
agreed well with the theory. However, the available 
data were insufficient to examine critically the 
effects of two important properties on k& — the 
thermal conductivity of the fluid and the pressure. 
The objectives of the present study are twofold: 
(1) to measure the effect on k% of (a) the pressure 
for gases, and (b) the thermal conductivity of the 
fluid k, for both gases and liquids; and (2) to deter- 
mine for rocks of different porosities and permea- 
bilities the apparent solid conductivity k, and the 
consolidation parameter b,D,/k , required in apply- 
ing the theoretical equations. 

Measurements extrapolated to zero pressure are 
helpful in evaluating k, and the consolidation pa- 
rameter. However, the extrapolations must be based 
upon data obtained at low pressures because in 
smal! pores the molecular conductivity of the fluid 
may be unusually depressed. This occurs when the 
mean free path of the molecules is of the same 
order of magnitude as the pore diameter, that is, in 
small-diameter pores at low pressures. The phenom- 
enon has been studied recently by Schotte® in beds 
of unconsolidated fine particles. With particles in 
the size range of 45 to 200 microns, a significant 
reduction of kg with pressure was found even at 
absolute pressures greater than 1 atm. In the sand- 
stones employed in the present investigation, pore 
diameters lower than those experienced by Schotte 
are likely, so that the effect would presumably be 
greater. 


SCOPE OF MEASUREMENTS 


Stagnant conductivities were measured for four 
types of sandstones using a total of seven samples. 
The characteristics of the samples, including po- 
rosity and permeability values, are given in Table 
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i. 
Four gases were employed, data being obtained 


with nitrogen, helium and carbon dioxide at pres- 
sures above atmospheric; all sub-atmospheric data 
were determined using air. To provide a reasonably 
wide range of fluid thermal conductivity, measure- 
ments also were made using four fluids — normal 
heptane, methyl alcohol, a solution of 79.8 weight 
per cent glycerol and 20.2 weight per cent water, 
and pure water. The range of kg values obtained 
with these seven fluids is approximately 0.01 to 
0.36 Btu/(hr ft °F). 

The mean temperatures existing in the rock 
samples during the runs are given in Table 2 for 
each fluid. 


EXPERIMENTAL 


The apparatus for measuring the stagnant con- 
ductivity is shown in Fig. 1. The test section con- 
sisted of two cylindrical cores, each | in. in diam- 
eter and about 1 3/8 in. in length. To reduce longi- 
tudinal heat transfer, one additional core sample 
of the same dimensions was used on each end of 
the test section. The conductivity was determined 
in the radial direction of the cylindrical cores, that 
is, in the direction parallel to the plane of the 
sandstone formation. 

To introduce a radial temperature gradient, an 
electric heater was made by winding 30-gauge 
nichrome wire around a flexible rod. This rod was 
inserted into holes 0.093 in. in diameter drilled 
through the central axis of the core samples. The 
power input in the test section was determined 





TABLE 1 — POROUS ROCK SAMPLES USED 


Boise Sandstone from Idaho Outcrop 
Yellowish gray, medium-grained, argillaceous sandstone; 
micaceous, feldspathic, ‘‘pepper and salt’’ appearance; mas- 
sive, 
BO-1: porosity 0.248, permeability 590 md 
BO-2: porosity 0.256 

Bartlesville Sandstone from an Oklahoma Outcrop 
Moderate yellowish brown, fine-grained, argillaceous sand- 
stone; minor muscovite flakes, altered iron minerals, massive. 
BA-2: porosity 0.203, permeability 28 md 
BA-5: porosity 0.220 

Berea Sandstone from Ohio Outcrop 
Very light gray, fine-grained, clean, well sorted sandstone; 
massive. 
BE-1: porosity 0.189, permeability 170 md 
BE-2: porosity 0.185 

Rangely (Colorado) Reservoir Rock 
Olive gray, fine-grained, clean, well compacted sandstone. 
RA: porosity 0.117, permeability 18 md 

















TABLE 2 — RANGE OF MEAN TEMPERATURES OF SAND- 
STONE SAMPLES DURING RUNS 
°F 
Carbon Dioxide. 
Nitrogen. . 140 to 190 
Helium . . ° 90 to 100 
Normal Heptane . oe ee ee - 110 to 140 
Mey Aicene) « «6 6 kt ll Uw 90 to 110 
Giycersl®. « « «© » © © © «= «© © © » « MeOter too 
Se a ee ee ee ee ee 
* 79.8 wt. per cent glycerol and 
20.2 wt. per cent water 
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140 to 200 


with a Weston Wattmeter (Model #432) with an 
accuracy of 0.5 w. The measurement was made by 
attaching electrical leads to the heater wire across 
the length of the test section (approximately 2 3/4 
in.) Fine sand was used between the test cores and 
the inner wall surface to adjust the position of the 
cores, as illustrated in Fig. 1. Cooling water was 
passed through a jacket surrounding the sample 
when runs were made with gases. 

Copper-constantan thermocouples (30 B. and S. 
gauge) were used to measure the radial temperature 
gtadient. The thermocouples were inserted in holes 
drilled in the end surface of the separate core 
sample as illustrated in Fig. 1. The holes, 0.027 
in. in diameter, were drilled to a depth of 0.17 in. 
Three or four thermocouples were used across each 
of several radii, and their radial location was 
measured carefully to the nearest 0.1 mm. 

For runs made with gases as the stagnant fluid, 
the temperatures were measured using a recorder 
with an expected accuracy of 0.5°F. The fluid was 
alternately added and evacuated from the test sec- 
tion several times before making thermal measure- 
ments. Fluid pressures on the sample at atmospheric 
pressure and above were measured with a mercury 
manometer or a Heisse gauge. Sub-atmospheric pres- 
sures were determined with a vacuum gauge. 

A somewhat different procedure was employed 
with liquids. The solid cores, dried completely 
with an infra-red lamp, were first placed in a con- 
tainer half-filled with the required liquid. The re- 
maining gas in the pores of the core was removed 
with a vacuum pump. Upon removal from the con- 
tainer, the samples were covered with scotch tape 
to prevent evaporation of fluid. Then, as quickly 
as possible, the thermocouples were inserted. The 
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FIG. 1 — EXPERIMENTAL APPARATUS FOR GASES. 
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test element was then placed in a container filled 
with the required liquid. The entire unit was in- 
serted in an ice-water, constant-temperature bath. 
Temperatures were determined for the runs with 
liquids using a precision portable potentiometer. 


EVALUATION OF EXPERIMENTAL PROCEDURE 


If longitudinal heat transfer in the core samples 
is negligible, the stagnant conductivity can be 
calculated from the measured radial temperature 
profile by the expression, 

».. _. £iae 

Re = IaL dT je @ 3 owe SE 

Typical measurements for several fluids in Core 
Sample BE-1 are shown in Fig. 2. Within the pre- 
cision of the data, a straight line relationship is 
observed. This suggests, according to Eq. 1, that 
the variation of k% with temperature across the 
radius of the sample is not sufficient to be measured 
by the experimental procedure used. Hence, a single 
mean value for k% is calculated from Eq. 1 using 
the slope of the line obtained from plots similar to 
those illustrated in Fig. 2. 

To test the assumption of negligible longitudinal 
heat transfer, several preliminary measurements 
were made using Sample BA-5 with air. First, radial 
temperature profiles were observed and compared 
at the three contact surfaces (the ends of the cylin- 
drical core samples) of the four core samples mounted 
in the apparatus (Fig. 1). These results indicated 
that approximately 0.4 per cent of the total energy 
input from the heater flowed out in the longitudinal 
direction. 

As the energy input from the heater decreases, a 
point is ultimately reached when longitudinal heat 
losses will be significant. Preliminary measure- 
ments indicated that this energy level was about 
10 w. All the runs were made at values of Q greater 
than this limiting value. . 
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FIG. 2 — EXAMPLES-TEMPERATURE DISTRIBUTION 
IN CORE SAMPLE, 
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To test further the significance of longitudinal 
heat transfer, measurements were carried out when 
the ends of the apparatus were maintained at differ- 
ent temperatures with nichrome wire heaters. The 
energy input to the central heater was maintained at 
15 w, leading to a mean core temperature of about 
150°F. Under these conditions, varying the temper- 
ature of the ends of the apparatus from 70° to 200° 
F did not affect the stagnant conductivity. 

From these preliminary tests it appeared that the 
equipment was satisfactory for measuring reasonably 
accurate values of k%. Considering errors in energy 
input, temperature and position measurements, and 
errors in applying Eq. 1, it is believed that the max- 
imum uncertainty was less than 10 per cent. 


RESULTS 


EFFECT OF PRESSURE ON k§ 

The data showing the effect of pressure on k@ for 
gases are plotted in Figs. 5 and 6 for pressures 
greater than 1 atm. Fig. 5 demonstrates that different 
core samples from the same sandstone formation 
(Bartlesville sandstone) do not exhibit significant 
differences in stagnant conductivity. On the other 
hand, the thermal conductivity of the gas is a sig- 
nificant factor; i.e., the results for helium are ap- 
preciably higher than those for nitrogen and carbon 
dioxide. In Fig. 6, data are plotted for three different 
types of sandstones, and here large differences are 
observed. 

The small increase in k% with pressure at high 
pressures, shown in Figs. 3 and 4, is probably due 
to the small increase in the normal thermal conduc- 
tivity of the gas with pressure. However, the more 
rapid change in kQ with pressure at low values of 
P must be the result of another factor because the 
normal thermal conductivity is essentially constant 
at these pressures. This effect is more pronounced 
at sub-atmospheric pressures as demonstrated by 
the curves in Fig. 5 for air. The decrease in k% 
with pressure is most pronounced for the sample 
with the lowest porosity, Rangely sandstone (RA). 
This suggests that the fluid conductivity is being 
reduced by the small-diameter pores in the core 
samples. Schotte® has adapted equations based 
upon kinetic theory to show the effect of pressure 
on the apparent value of k, in packed beds of fine 
particles (45 to 200 microns). One of his equations 
can be reduced to the following form to illustrate 
the effect of pressure P and pore size 6. 

ke 

"s-TeC(TAb cD 
Here, k*, is the normal thermal conductivity of the 
gas. For a given core sample at a fixed temperature, 
Eq. 2 indicates that k, approaches the normal 
conductivity as P approaches a large value and 
that k, approaches zero as P tends toward zero. 
This suggests that the stagnant thermal conductiv- 
ity of the core (k%) would approach a constant 
value, equal to the contribution of the solid phase, 
as the pressure approaches zero. The S-shaped 
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behavior of the curves in Fig. 5 is in agreement 
with this reasoning. The significant conclusion 
here is that the solid contribution to the stagnant 
conductivity (k%), cannot be obtained by extrapol- 
ating to P =0 data at above atmospheric pressure, 
such as shown in Figs. 3 and 4. Rather, measure- 
ments must be made under vacuum conditions to 
obtain (k°,), by extrapolation. The values of this 
parameter, as determined from Fig. 5 for the four 
types of sandstone, are given in Table 3. 


EVALUATION OF CONSOLIDATION PARAMETER 
In Ref. 4, the following equation was proposed to 
predict the stagnant conductivity of porous rocks 


Porosity Permeability CO2 No 
BO-! 0248 590 ad a ° 
BE-! 0.189 170 ama A my 


RA 0.117 18 ad ra ° 
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The objective here is to use this expression and 
the experimental values of (k@), in Table 3 to 
determine the relationship between the apparent 
thermal conductivity of the solid phase k, and the 
consolidation parameter hyD,/k,. Both are charac- 
teristic properties of the sandstone sample. This 
is accomplished by considering Eq. 3 at zero pres- 
sure where k, approaches zero. From the equations 
in Ref. 4, the following limiting results are readily 
obtained. 


Lim (1 ke\ _ 

kg +0 (: it) =< a oe oe oe oe . (4) 
and 

imam, 4h © 4 et wom wie, @ & 6 pe) 
kg +0 


The value of «, the void fraction of the original 
bed of unconsolidated particles prior to consolida- 
tion, is suggested in Ref. 4 to be 0.476. This is 
the maximum void fraction in a uniform arrangement 
of spherical particles in contact. Using this result 
and Eqs. 4 and 5, the general expression(Eq. 3) for 
k°./k reduces at zero pressure to the form, 


0B) 


This equation provides a relationship between k, 
and Dphy/ks. If ks were known, D, hy/k, could be 
determined from (k%)> and Eq. 6. However, kg is 
not necessarily. the thermal conductivity of the 
pure solid phase (for example, quartz for sandstone 
cores). It may differ from this due to the non-uniform 
geometry of the solid phase and due to impurities 
between the grains of solid particles and the ce- 
menting material holding the grains together. 

This problem is avoided by utilizing the experi- 
mental data for the effect of kz on k%o, along with 
Eq. 6, to ascertain separate values for k, and 
Dyhp/ks- The experimental data for the several 
fluids and sandstone samples are shown in Figs. 6 
through 8. The curves on the figures represent Eq. 
3 plotted for certain values of k, and Dphy/ks. 
These specific values were obtained in the follow- 
ing way. 

1, Sets of consistent numbers for k, and Dy hy/ks 
were first determined from Eq. 6. 

2. The most appropriate set was then obtained 
for each type of sandstone by comparison with the 
experimental data in Figs. 6 thtough 8. 

Figs. 6 and 7 include calculated curves for three 
sets of k, and hyD,/k,, designated on the curves 
by the numerical value for k,. These results indi- 
cate that it is not possible to determine k, more 
closely than approximately + 0.2 Btu/(hr ft °F). 





(6) 





TABLE 3 — THERMAL PROPERTIES OF SANDSTONE 





Sandstone (k Qo ks Pphp 
Sample Btu/ft he °F Btu/ft hr °F ks 
BO 0.88 1.8 0.965 
BA 0.69 3.2 0.159 
BE 0.77 4.6 0.103 
1.22 4.6 0.105 
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The values of the two parameters so evaluated are 
shown in the last two columns of Table 3. 

It is interesting to note that these parameters 
are essentially the same for the BE and RA sand- 
stones. The difference in stagnant conductivity 
between the two types of rock must be due primarily 
to the difference in porosities. Visual examination 
of the Boise (BO) sandstone suggests that it has a 
markedly different structure from the other types. 
It seems to be well consolidated with somewhat 
larger pore diameters than those associated with 
the other sandstones. This micaceous and feld- 
spathic structure is probably responsible for the 
low value of k, (Table 3) determined for the Boise 
sample. 


CORRELATION AND PREDICTION 
OF STAGNANT CONDUCTIVITIES 


For Eq. 3 to be a useful prediction method for 
k°e, it is necessary that values of k, and Dphy/ks 
be available. At present these values cannot be 
obtained except by experimental means. Until data 
such as that shown in Table 3 are available for a 
large variety of rock samples, the method is not 
valuable for prediction purposes. However, the 
agreement in Figs. 6 through 8 between experimental 
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data and Eq. 3 suggests that the procedure correctly 
accounts for the effects of properties of the rock 
and fluid, such as void fraction and fluid conduc- 
tivity kz. Hence Eq. 3 is useful now as a means of 
predicting the effect of these properties on ke and 
requires no experimental data for this use. The 
method is also valuable for correlating experimental 
ke information. 

Of the data in the literature, only that of Zierfuss 
and Vliet® included information on the effect of 
different fluid conductivities. Their conductivity 
results for air and liquid water showed increases 
in k%& with kg of the same magnitude as indicated 
in Figs. 6 through 8. Data were presented for only 
these two fluids with known thermal conductivities, 
and no information was given at low pressures. 
Hence, it was not possible to evaluate (k%), and 
make a quantitative analysis of their results accord- 
ing to Eq. 3. 


NOMENCLATURE 
b = distance over which conduction takes 
place (i.e., pore diameter), ft 
C = constant, (psi) (ft)(1/°R) 
Dphyp/ks = dimensionless consolidation parameter 


for expressing the heat transfer in the 
solid phase (see Ref. 4) 
i stagnant conductivity, Btu/(ft)(hr) (°F) 


asymptotic value of stagnant conductiv- 
ity corresponding to zero thermal con- 
ductivity of gas phase (i.e., zero pres- 
sure), Btu/(ft) (hr) (°F) 
kg = apparent thermal conductivity of gas, 
Btu/(ft) (hr) (°F) 
k*, = normal thermal conductivity of gas, Btu/ 
(ft) (hr) (°F) 
k, = apparent thermal conductivity of solid 
phase, Btu/(ft) (hr) (°F) 
L = length of test section, ft (approximately 
2 3/4 in.) 
= absolute pressure, lb/sq ft 
= radial distance, measured from center 
of core, ft 
= energy input in the test section, Btu/hr 
absolute temperature, °R 
= ratio of effective thickness of the fluid 
film adjacent to the contact surface of 
two solid particles in the original 
packed bed (prediction of ¢ is de- 
scribed in Ref. 4) 
@ = void fraction of original packed bed of 
unconsolidated particles 
’ = porosity of sandstone 
dlinr _ 
dT 


a2 
I 


% 
U 


q oO 
Nl 


slope of temperature vs radial distance 
line plotted as Inr vs T (Fig. 2), °R7} 
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. ABSTRACT 


A rigorous study was made of the behavior of 
reservoirs composed of horizontal layers, uncon- 
nected except at the well and filled with a compres- 
sible fluid. The report is presented in two parts. 
Part I deals with the practical implications of the 
results obtained; Part II contains the mathematical 
derivation of the solution. The main portion of Part 
I is concerned with application of the results to 
reservoirs composed of two layers. Relative rates 
of depletion of the layers are studied, and it is 
shown that differential depletion between layers 
exists in the transient stage of the reservoir, which 
generally is of an order of magnitude longer than 
the transient stage in a single-layer reservoir. In 
the transient stage, the more permeable layer is 
depleted faster than the less permeable. However, 
the reservoir approaches a steady-state situation 
where both layers contribute equally to production. 
Theoretical build-up curves from such reservoirs 
are presented, and the influence of a skin effect 
and flow into casing and tubing are also evaluated. 
The build-up curves are predicted to have, first, 
the familiar logarithmic straight-line section and, 
subsequently, a rising and flattened section. The 
results found are compared with an extension of an 
approximate theory, and it is found that this simpli- 
fied theory is applicable over most of the life of 
the reservoir. Extensions to reservoirs containing 
more than two layers are indicated, and examples 
of the performance and build-up curves for three- 
layered reservoirs are presented. The theoretical 
results on build-up curves were applied to field 
examples of such reservoirs, and satisfactory re- 
sults were obtained. From the build-up curve, the 
permeability-thickness product, the wellbore dam- 
age and the static pressure are obtained. However, 
it apparently is not possible to determine the prop- 
erties of the individual layers from the combined 
build-up curve. 
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PART I — PRACTICAL APPLICATIONS 


INTRODUCTION 


A large number of reservoirs have been found 
where the producing formation is composed of two 
or more layers of differing physical characteristics, 
such as permeability and porosity, and where also 
the thicknesses of the layers differ. Since in general 
a differential depletion should exist between the 
layers of such reservoirs, it would be advantageous 
to know the amount of the differential depletion 
under a given production history. It also would be 
valuable to know the pressure behavior, because it 
could be undesirable to produce a reservoir in such 
a manner that the pressure in one of the layers 
would have fallen below the bubble point of the 
oil at a time when another layer might still be 
virtually undepleted. 

The determination of static pressures in shut-in 
wells with slow build-up is another feature of im- 
portance in layered reservoirs. If the behavior of 
the reservoir differs appreciably from the predicted 
behavior of a single-layer reservoir, the type of 
analysis for build-up curves suggested for single- 
layer reservoirs is not applicable. 

Several authors have treated the problem of find- 
ing the pressure and production characteristics of 
multilayer reservoirs. Horner! treated the problem 
of n layers in an infinite field, the layers being 
connected only at the well. Tempelaar-Lietz? has 
presented an approximate treatment of the charac- 
teristics of a bounded reservoir composed of two 
layers of different permeabilities and equal thick- 
nesses. Horner gave no information on differential 
depletion between layers, and the treatment of 
Tempelaar-Lietz was not applicable for prediction 
of pressure build-up behavior. In addition, it was 
thought that a more rigorous treatment than the 
one by Tempelaar-Lietz would be valuable in de- 
termining the range of validity of this solution. 
The present study, therefore, was made of a 
bounded reservoir composed of two or more layers, 
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these layers being unconnected except at the 
well, 


BASIC ASSUMPTIONS 


The reservoir to be studied is shown in Fig. 1. 
Each layer of the reservoir is assumed to be homo- 
geneous and isotropic, and filled with a fluid of 
small and constant compressibility c. The permea- 
bility k, oil-filled porosity @ and thickness h of 
each layer are designated by the corresponding 
symbol with the subscript 1, 2, ...., m, depending 
on which layer is referenced. The viscosity of the 
oil is assumed to be constant and is denoted by p. 
The reservoir is initially at a uniform pressure pp; 
and at all times ¢t > 0, the fluid is produced in such 
a manner that the production rate q7, measured at 
initial reservoir conditions, is held constant. 

Under these assumptions, the pressure at any 
point must satisfy the well known equation for flow 
of a fluid of small but constant compressibility 
given in Part II. The solution of this equation for 
the proper boundary conditions is obtained with the 
aid of the Laplace transform. 


DISCUSSION OF RESULTS 


FLOWING-WELL PERFORMANCE 
For all but very small times, the solution to the 
differential equation leads to the following results: 


bj — Pw(t) = 
ar|__ = "w + Y(t) 
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are directly proportional to qr. 
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FIG. 1 — THE RESERVOIR. 


For the case of two layers, these equations can 
be rewritten as 
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For early times at which the influence of the 
boundary has not been felt, the expression 


bi Py t) 











qrp/4mkb 
kyhy ln Pr Her y + koh In Pau = 
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gives the correct pressure behavior at the well, y 
denoting Euler’s constant. The lower limit of the 
time range of validity of this equation is small 
enough to include most times of practical interest. 

From Eqs. 3 and 4 it can be seen that a steady 
state is finally attained in the reservoir, during 
which the rate of change of the well pressure with 
respect to time (i.e., dp,/dt) is constant. Further- 
more, at this time the production rate from each 
layer becomes constant, and it then follows that 
the pressure everywhere within the reservoir varies 
linearly with time. 

Solutions for two sample cases for two sets of 
parameters are presented in Figs. 2 and 3. The 
curves labeled I correspond to the parameters 


ki P1 
ae ial 
ko $2 
a “€- 2,000 ; 
bo Tw 
whereas, the curves labeled II correspond to the 
parameters 
a aR ; or =1 ’ 
k2 $2 
a ' “e - 2,000 ; 
ho Tw 


Fig. 2 shows the pressure behavior at the well. 
Eq. 5 was used in calculating this behavior for tp 
< 10°. During this first period, the pressure de- 
creases logarithmically with time. Physically, this 
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time range corresponds to the time during which 
the well has not detectably influenced the pressure 
at the outer drainage boundary. Hence, the pressure 
decline at the well is identical to the decline of a 
well draining an infinite reservoir. The time at 
which the influence of the drainage boundary is 
first detected at the well is characterized by the 
break of the pressure curve away from the straight- 
line portion of the curve. Now the reservoir pro- 
duces as if it were composed of one bounded layer 
and one infinite layer because the influence of the 
: boundary has been reflected only in the more per- 
) meable layer. After some time, the less permeable 
layer also feels its boundary. From then on, the 
F reservoir behaves as a fully bounded reservoir. 
Steady state, where the relation between p,, and t 
) would be a straight line if plotted on a linear scale, 
is reached some time thereafter. The dotted lines 
in Fig. 2 are the approximations given by the method 
of Tempelaar-Lietz and will be discussed in a later 
section. : 

The magnitude of the time necessary to reach 
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pared with the corresponding time for a single-layer 
reservoir, Matthews and Brons? have shown that for 
a cylindrical reservoir steady state is reached at a 
dimensionless time 


kt 

pper2 

If r./r,, = 2,000, then steady state is reached at 
kt 

pycri, 


In a two-layer reservoir, the dimensionless time 
needed to reach steady state varies with the prop- 
erties of the two layers; the average of the cases 
studied in this report is about 





= 0.3 , 





= 12x 10° 


kt 
ducr2 = 8 x 10” 


Hence, the time needed to reach steady state is of 
the order of 50 times as great as the corresponding 
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FIG. 3— FRACTIONAL PRODUCTION RATES FROM THE MORE PERMEABLE LAYER OF 


TWO 2-LAYER RESER- 








this lies mainly in the changing rates of production 
from the layers in the multilayer case. 

The behavior of the production rate for these two 
cases is shown in Fig. 3, in each instance for the 
more permeable layer. The ordinate in the graph is 
qi(t)/qr, the fractional production rate from this 
layer. The dotted lines again are the approximations 
of Tempelaar-Lietz. The fractional production rate 
from the other layer qo(t)/q7 can be obtained from 
the relation 


q(t) , q(t) 
qT qT 


At early times, the production rate from the more 
permeable layer is increasing; it approaches asymp- 
totically the value q4(t)/q7 = kyby/(kyby + kob2), 
which behavior corresponds to the case of a well 
draining an infinite reservoir. After the initial in- 
crease, the production rate begins to decrease, a 
phenomenon which corresponds to the physical fact 
that the influence of the boundary in the more per- 
meable layer has been reflected on the behavior of 
the reservoir. At some later time, the effect of the 
boundary is also felt in the less permeable layer; 
subsequent to this time, the curves are seen to 
approach asymptotically a constant value, this value 
being 93(t)/q7 =b1¢1(h1 1 + 5262). This value is 
the fractional production rate of Layer 1 at steady 
state. 

The ratio by ¢4/(by¢1 + h2¢2) is equal to the 
ratio of the oil-filled volume of Layer 1 to the oil- 
filled volume of the entire reservoir; hence, at 
steady state, each layer produces at a fractional 
rate equal to the fraction of oil-filled volume it 
contains. 

With regard to fractional depletion between the 
layers, if each layer always produced at a fractional 
rate proportional to the oil-filled volume it contains, 
there would be no differential depletion. Hence, 
during production of a reservoir, differential deple- 
tion takes place previous to the time when steady 
state is reached. 

The total production rate q7 from the reservoir 
is an important factor in the differential depletion 
of the reservoir. The time required to reach steady 
state is essentially independent of gz, but the 
cumulative production from each layer is directly 
proportional to qr. If Q(t) and Qo(t) denote the 
cumulative production at the time ¢ from Layers 1 
and 2, respectively, 





t 
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Also, the differential depletion AD(t) per unit thick- 
ness is then defined as 


Q(t) Qylt) 


AD(t) = by 








But we can also write 
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Since the ratios q,(t)/q7 and q2(t)/gy are inde- 
pendent of q7, the differential depletion is directly 
proportional to gz. For example, if the maximum 
differential depletion in a reservoir produced at a 
rate q7 is AD and if the reservoir had been pro- 
duced at a rate 2 qy, the maximum differential de- 
pletion would have been 2AD. 

The results from another sample case are pre- 
sented in Figs. 4 and 5. The properties of the 
layers are ky/kz = 4, by/hg =0.05 and 44/2 = 2. 
The results are presented for four different values 
of re/ty —Te/ty =-2,000, 1,000, 100 and 10. Fig. 
4 shows the well pressure behavior as a function 
of production time; the straight-line portion again 
corresponds to the time range when the influence 
of the boundary has not been felt at the well, and 
the deviation from this straight line corresponds 
to the influence of the boundary. 

The times at which the pressure-decline curves 
deviate from the straight-line portion increase with 
with the dimensions of the reservoir. This is be- 
cause the time required (t,) for the pressure dis- 
turbance caused by the opening of the well to reach 
a certain point (r) away from the well increases 
with the distance of the point from the well. A 
good approximation is given by t, « r2, 

Fig. 5 shows the fractional production rate from 
the more permeable layer for several ratios re/ry, 
for the properties just mentioned. 


THEORETICAL PRESSURE BUILD-UP CURVES 

Since the behavior of the reservoir is described 
by linear differential equations, the behavior during 
build-up can be obtained by superposition. Fig. 6 
shows a typical theoretical pressure build-up curve 
obtained from a two-layer reservoir. As in a single- 
layer reservoir, there is an initial straight-line 
section, AB. After the straight-line portion, the 
build-up curve levels off (BC), this leveling off 
corresponding in a single-layer reservoir to the 
pressure’s having almost reached its average pres- 
sure. However, in a multilayer reservoir the pressure 
again rises (CD) and then finally levels off at the 
average pressure (DE). The rise in CD is due to 
the repressuring of the more depleted, permeable 
layer by the less depleted, less permeable layer. 

Several other examples of pressure build-up curves 
are presented in Figs. 7 through 10, the parameters 
for each reservoir being indicated on the same page 
as the graph. Figs. 7 and 8 present pressure build- 
up curves, plotted both against At, the shut-in time 
and against At/(t + At). Figs. 7 and 8 are ata 
closed-in time when the reservoir is still in the 
stage where the influence of the boundary has not 
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FIG. 5—FRACTIONAL PRODUCTION RATES FROM MORE PERMEABLE LAYER OF A TWO-LAYER RESERVOIR 


FOR DIFFERENT RATIOS OF r,/r,,. 


been felt during production. Fig. 9 is a build-up 
curve for the same reservoir at a closed-in time 
when the influence of the boundary has been felt. 
From these curves, it can be seen that the entire 
portion CDE (as in Fig. 6) expands greatly in ver- 
tical scale with increasing production time. 


The magnitude of the pressure rise during the 
portion CDE of the pressure build-up curve depends 
on the contrast of the properties of the layers. If 
the two layers are nearly equal in permeability, 
the pressure rise in this portion will be small (see 
Fig. 10), whereas if the two layers differ widely in 
properties, the pressure rise will be considerable 
(see Fig. 9). 

For equivalent information, a well in a multilayer 
reservoir must remain shut in longer than a well in 
a single-layer reservoir, the ratio of required shut- 
in times being roughly equal to the ratio of the 
times needed to attain steady state in the two res- 
ervoirs, It was seen that this latter ratio is large 
(of the order of 50); it then follows that wells pro- 
ducing from multilayer reservoirs may have to re- 
main shut in for considerable times to provide use- 
ful pressure build-up data. 

A reason for the difference in type of pressure 
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build-up curves obtained from two-layer and single- 
layer reservoirs is that in a two-layer reservoir 
flow occurs from one layer to the other during build- 
up. There is differential depletion between the 
layers, the more permeable layer having been de- 
pleted more than the less permeable layer; when 
the well is shut in, there must be flow from the 








At/(t+ At) 


FIG. 6 — THEORETICAL PRESSURE BUILD- 
UP CURVE FOR TWO-LAYER RES- 
ERVOIR, 














less permeable layer to the more permeable layer 
for the pressure to stabilize over the entire reser- 
voir. 


SKIN EFFECT 
In the study of skin effects, the skin factor S; in 
each layer was defined by means of the equations 
ee J ee © ne Gin 


where ps; denotes the formation pressure in the jth 
layer. 


The mathematical analysis of the problem showed 
that the effect of the skin could be represented by 
a fictitious well radius in each layer. Eqs. 1 and 2 
were modified to include these characteristics. 

The effect of introducing a fictitious well radius 
can best be illustrated for the case of a single-layer 
reservoir producing at a constant rate. As shown in 
Fig. 11, the skin will cause a discontinuity in the 
pressure at the well. The behavior of the pressure 
in the vicinity of the well is logarithmic; if the 
curve is extrapolated logarithmically (the dotted 
line in Fig. 11), there exists a point r‘%, such that 
p(r*) = Py. Then r%, can be considered the ficti- 
tious well radius representing the skin. The rela- 
tion between r,, and r., is given by 


wate ok, Pe es. 


Several curves showing the influence of skin 
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FIG. 7— PRESSURE BUILD-UP CURVE WITH At/(t + At) 
AS ABSCISSA. 
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FIG. 8—PRESSURE BUILD-UP CURVE WITH At AS 
ABSCISSA. 





effects are presented in Figs. 12 and 13. The char- 
acteristics of the reservoir studied are 
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In each graph, Curve a refers to the case where 
no skin was present in either layer (i.e., 5; = 0 and 
Sq = 0); Curve b refers to the case where the skin 
effects were taken to be S; = 10 and Sz = 1; and 
Curve c refers to the case where the skin effects 
were taken to be Sj = 1 and Sy = 10. 

A comparison shows that for Case b there was 
less differential depletion than for Case a, where 
the reservoir had no skin, and that for Case c there 
was more differential depletion than for Case a. 
These facts are readily explained; if a high resist- 
ance to flow is placed in the more permeable layer, 
it will tend to equalize the flow from the layers 
and, thus, will reduce differential depletion. 

Pressure build-up curves for these cases and for 
cases with no skin have the same sort of shape. 
However, the effect of the skin on pressure build- 
up curves is different for multilayer and single- 
layer reservoirs. For a reservoir composed of a 
single layer, the effect vanishes after a short in- 
terval of shut-in time; on the other hand, in a mul- 
tilayer reservoir the effect of the skin is present 
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during the entire range of shut-in time. This is true 
because, at the time the well is closed in, differ- 
ential depletion exists between the layers; during 
shut-in, flow must then occur from one layer to 
another through the wellbore, and the skin acts as 
a resistance to that flow. Hence, the influence of 
the skin is felt during the entire pressure build-up. 


FLOW INTO CASING AND 
TUBING DURING BUILD-UP 


To evaluate this effect, it was assumed that the 
flow from the formation did not drop discontinuously 
to zero at the instant the well was closed in, but 
that the flow decreased in an exponential fashion; 
that is, 


q-( At) = gre" tt , ak eB cene le (D) 


as discussed by van Everdingen* for the single- 
layer case. 

Results show that after-production deforms the 
straight-line portion of the pressure build-up curve 
in the same fashion that it influences the curve for 
a single-layer reservoir. Fig. 14 shows a build-up 
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FIG. 11—FICTITIOUS WELL RADIUS TO 
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FIG, 12—PRESSURE-DECLINE CURVES SHOWING IN- 
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curve both with and without after-flow. The two 
curves overlap for all shut-in times for which the 
flow due to after-production is essentially zero. 


EXTENSION TO THREE-LAYER RESERVOIRS 

Several examples of three-layer reservoirs were 
studied; results show that, qualitatively, their 
behavior is like that of a two-layer reservoir. Dif- 
ferential depletion exists between layers, the dif- 
ferential depletion being directly proportional to 
the total production rate. The shape of the pressure 
build-up curves is the same as that for two-layer 
reservoirs. Figs. 15 and 16 give examples of re- 
sults obtained for three-layer reservoirs. 


COMPARISON OF RESULTS WITH 
THE TEMPELAAR-LIETZ THEORY 


The previously mentioned work of Tempelaar- 
Lietz presents an approximate theory of the behavior 
of a two-layer reservoir where the layers had dif- 
ferent permeabilities, equal thicknesses and equal 
porosities. The simplifying assumption made by 
Tempelaar-Lietz in his treatment was that the pro- 
duction rate from each layer was proportional to 
the difference between average pressure in that 
layer and the well pressure, the constant of propor- 
tionality being a function of the parameters asso- 
ciated with that layer. 

The theory presented by Tempelaar-Lietz easily 
can be extended to include differences between the 
layers in thicknesses and porosities, as well as in 
well radii. The formulas obtained from this treat- 
ment are given at the end of this section. 

The results obtained by the Tempelaar- Lietz 
theory were compared with the exact solutions 
obtained in this study. It was found that the Tem- 
pelaar-Lietz solutions are very good approximations 
to the true solutions for pressure decline and rate 
behavior for all but early production times. To be 
more exact, the Tempelaar-Lietz solutions do not 
express the true behavior for early times (when the 
reservoir can be considered to behave like an in- 
finite reservoir) and for times shortly thereafter. 
After the ‘influence of the boundary is felt, the 
Tempelaar-Lietz solutions become applicable with 
good approximation. A comparison of results is 
presented in Figs. 2 and 3, where both the true and 
approximate solutions are presented for one example. 
The approximate solutions are labeled with the 
subscript a. 

The Tempelaar-Lietz pressure-decline equation 
cannot be used as an approximation to the pressure 
build-up curves because the pressure values are 
needed at times at which the Tempelaar-Lietz so- 
lution is inapplicable. However, the portion CDE 
in Fig. 6 can be described by the Tempelaar-Lietz 
equation, and the method for finding average pres- 
sures described later will make use of this fact. 

The generalization of the Tempelaar-Lietz theory 
to the case of unequal properties of the layers gives 
the following formulas. 
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and q(t) can be obtained from q2(t) = qT — 91(t). 
In these equations, A; is the recovery per unit 








; pressure drop of Layer j and i; is the productivity 
index of Layer j, as defined by 
{ Aj = cpyr2by 5.2... CA) 
: Qn keh; 
ee 
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APPLICATIONS 


RELATIVE DEPLETION BETWEEN LAYERS 
One application of the results of this study lies 
a in the prediction of relative depletion between 
. layers as a function of withdrawal rate. 
If the rate of withdrawal q is not constant, the 
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solutions for varying rate can be built up by appro- 
priate superposition. The following example should 
make this clear. 
Example 

Let a well which penetrates two layers be pro- 
duced at a rate of 100 B/D for 30 days, shut in one 
day, then produced at a rate of 50 B/D. Find the 
differential depletion between layers at the end of 
60 days. The physical constants of the reservoir 
are as follow: ky = 0.05 darcy, kz = 0.005 darcy, 
hy = 25 ft, by = 75 ft, c = 10-5 psi-!, p =0.5 cp, d 
= 0.2, re = 500 ft and r,, = 1/4 ft; t’ denotes pro- 
duction time in days. 


Then 
2 
Ay = sorert = 6.99 bbl/psi, 
A2 = 20.97 bbl/psi, 
F _ _2mkyhy 1 1 
x “2 toa 
B Tw 4 

= 2.583 B/D x psi , 

i2 = 0.775 B/D x psi, 
and 


pi -Pi = a6 + 0.652 q(1-e-0.1137 t’) , 
For the production schedule outlined, we have at 
the end of 60 days, 


p; -?, = ie +0.652(100) [1 —e -0-1137(60)] 


is —— —0.652(100)[1—e-0- 1137(60-30) ] 
.96 


+ 2KG0=3) | 9.652(50)[1 -e-0-1137(60-31)] 

27.96 
192 psi. 

—_— 100(30) + 50(29) _ 6.99 
20.97 20.97 
= 148 psi. 
OQ, =6.99(192) = 1,342 bbl or 53.6 bbl/ft sand, 
Qg = 20.97(148) = 3,103 bbl or 41.4 bbl/ft sand. 
For a rate of 300 B/D, the corresponding pressures 
at the same cumulative total production are 
bi — P1 = 319 psi, 
bi — p2 = 106 psi 
and the cumulative production per ft, 
O;/h, = 89.1 bbl/ft sand ’ 
Qo/hg = 29.6 bbl/ft sand . 

The effect of rate on the differential depletion 
thus is emphasized. These equations can also be 
used to calculate how near the average pressures 
in each layer approach one another during shut-in. 
The effect of a skin in each layer is taken into 
account by changing r, in each layer according to 
the manner of Eq. 8. 

PRESSURE BUILD-UP FOR MULTILAYER CASES 

As shown previously, the pressure build-up 
curves for multilayer cases have a “‘tail’’ which 
gives an essentially different shape to the build- 
up curve from that obtained for one-layer reservoirs. 
The method for obtaining the average pressure in 


tl 





(192) 


ll 


> 


51 





the drainage area of wells in such cases thus must 
be different from that for more homogeneous cases. 
The initial section of such a build-up curve is 
linear on a plot of pressure vs log [At/(t + At)] as 
shown in Fig. 6 where At is the shut-in time. Next 
a slight flattening may occur, then a rise and, after 
a long time, a final flattening. In the initial straight- 
line portion, a value for (ky, + k2h2) can be ob- 
tained from the slope by using 


qth/ 4a 
slope, psi 


kh) +kobo = x (14.7) (2.303) , 


where g7 is the total production rate measured at 
reservoir conditions at the time of shut-in. 

The applicable equation for the section of the 
build-up curve subsequent to the straight-line sec- 
tion, provided the well has been producing long 
enough before shut-in, is 


echt 


ote ee a 


where b and C are constants. Thus, a plot of log 
(p - p) vs At should give a straight line. One 
method of determining p is to assume a value for 
p, make this type of plot, and see whether or not 
a straight line is obtained. If not, a new value of 
p is assumed and the plotting is repeated. Fig. 17 
shows a series of these plots for different assumed 
values of p. The correct value is seen to be about 
2,750 psi because the lines on either side curve 
away from each other. 

It has been found that types of reservoirs other 
than multilayer reservoirs also show this same 
exponential behavior during build-up. Fig. 18 is a 
build-up curve from a fractured dolomite reservoir, 
the dotted line being the extrapolation to the aver- 
age pressure determined by the plot in Fig. 19. In 
Fig. 20 a build-up curve from a hydrafracked reser- 
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voir is presented, the average pressure having 
been determined with the aid of Fig. 21. Fig. 22 
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FIG, 21—SOUTH TEXAS FIELD (HYDRA- 
FRACKED RESERVOIR) OCT. 24 
TO OCT. 27, 1952. 





is a build-up curve from a multilayer* reservoir in 
California. The exponential-rise type of build-up 
appears applicable to all three types of reservoirs. 


CONCLUSIONS 


The present investigation has shown that the 
simple formulas developed by Tempelaar-Lietz and 
generalized in this report can be used successfully 
for predicting differential depletion between layers 
of a multilayer reservoir except at very early times. 
The simpler formulas can also be used to predict 
pressure decline at the well. A method for deter- 
mination of average pressures of these reservoirs 
is presented. This method can be applied to multi- 
layer reservoirs and appears to be applicable to 
reservoirs of other types which have the same type 
of build-up behavior. 


PART II — MATHEMATICAL DERIVATION 


ASSUMPTIONS AND BASIC EQUATIONS 


We assume that the reservoir to be studied is a 
bounded, horizontal cylindrical disk, with a com- 
pletely penetrating well at its center, that its 
thickness is small enough that the forces of gravity 
can be neglected and that it is bounded above, 
below and on the outside by impermeable bound- 
aries. The reservoir is divided horizontally into n 
parallel layers which are separated throughout the 
reservoir by impermeable boundaries through which 
no fluid can flow. The radius of the wellbore at 
the jth layer is denoted by ry,; and the exterior 
radius of the jth layer by re, ;, Each layer is assumed 
to be homogeneous, isotropic, and filled with a 
fluid of small and constant compressibility c; and 
constant viscosity p;(j = 1, 2,...., m). For each 





*For reservoirs composed of more than two layers, the appli- 
cable equation for build-up is of the type of Eq. 16, but it 
should contain n - 1 exponential terms if the reservoir contains 
n layers. Practically, however, it has been found that good 
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results generally can be obtained by using only one exponential, 
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FIG. 22—CALIFORNIA FIELD, PRESSURE BUILD-UP 
(MULTILAYER RESERVOIR) MAY, 1952. 


layer, k; denotes its permeability, ¢; its porosity 

and 4; its thickness. The reservoir is initially at 

a uniform pressure p;; and, at all times ¢ > 0, the 

fluid is withdrawn through the wellbore in such a 

manner that the total production rate g, measured 

at initial reservoir conditions, is held constant. 
For simplicity of notation we let 


kj 


PjMiCG 





ny = 





DS ae 
Vij 
ae a eres 


and define the pressure drops P;(r,t) by 

P,(1,t) = pj — pj(r,t) ju21,2,...,8 
where p,(7,t) denotes the pressure in the jth layer 
at a radial distance r and time t. Then, P; (r,t) must 
satisfy the equation 


Lo (, St). 4 % fo 
re a dae - Ey FeO ee ig Se 
« (€3) 


This system of equations must be solved subject 
to the following initial and boundary conditions. 
1, At ¢ = 0, the pressure drop in each layer is 
zero; i.€., 
P; (7,0) =0 (HT, 2iccy @ «+ « @® 
2. There is no flow across the outer boundaries, 
so that 


atr=re;,j=1, Bs eatig tt ora (3) 


3. The pressures in all layers at the well are 
equal and, hence, the pressure drops are also equal; 
i.e., 

P;(r,t) = P,,(t) atr= Tw,j» ] dd ee eee 
the function P,,(t) being independent of j. 
4, The total production rate from all the layers 


,n.. (A) 


53 








is constant and equal to 4g. 


Sai wae : B; r OP; es « 
i=1/ i=1'!? Or 


T=Tw,j 


where q(t) denotes the rate of flow from the jth 
layer into the well at time ¢. 


THE TRANSFORMED EQUATIONS 
AND THEIR SOLUTION 


Let X(z) = L[X(t)] denote the Laplace transform 
of X(t) where z is the complex variable x + iy. 
Applying this transform to both sides of the Eqs. 1 
and using Eq. 2, we obtain the system 


iz é au Zz 
r or Or nj 
. (6) 


which must be solved subject to the boundary 
conditions 


00 pwPhnccy 


wv! 


oP; 
—L=0 atr=ro,;, a a ae ee 
Or a 
» (7) 
P(r,z) = P,,(z), atr=rys, f=1,2,.+..,0 
(8) 
OP; , 
Ss ale — ee 
i=1 / Or * 2nz 


w,j 


Eq. 6 has general solutions of the form 


7 E: ‘z 
P(1,z) = A;Ky "nz + Bil, a 


and A; and B; are constants to be determined from 
Eqs. 7, 8 and 9. These conditions lead to the set 
of equations 


A; Ky (y; Vz) = Bly (y; V2) 
A;Ko(a;Vz) + B;Io(a;Vz) = P,,(z) .( 10) 


y B; a lAjK1(ajy2) — Bl (ajV/z)] = ea 


Letting 


Wo;(2) = Kola; V2)11(y; Vz) iv Ip(a; V2)K (yj; V2), 


Wy ;(2) 


we obtain by simple manipulation 


Ky (aj V2; Vz) - Iy(ajV2)Ky(y; Vz ), 





—" 1 
ra ™ as g PM ad ‘1 
j=1 7 Yo ,fz) 
and also 
P,(1,z) = 


Ko(-av2)h 1(%jV2) + Io ee ay/2 ol yjV2) _ 
Wo j(2) 7 





as Wz) . 
qj (2) = saa ae Pf os «+ § 5 
P07 


INVERSION OF THE SOLUTIONS 
In this section, we derive expressions for the 
inverses of the functions of Eqs. 11, 12 and 13 
suitable for numerical calculations at large times. 
Letting 


n Wy ,(z) 
F(z) = Vz > B;a; Ld caida 
j=1 Wo;(z) 
so that A , 
PAD ace ohn 
= 27z F(z) ’ 


we first note that V2; (z)/Uo (z) is a meromorphic 
function and that, therefore, the functions P,,(z), 
P;(r,z) and qj(z) are also meromorphic functions. 
From the stability of the physical system, it also 
can be deduced that these functions do not have 
poles with positive real part. 

Expanding 1/F(z) about z = 0, we obtain 








> Biy; (in i = 
2 1 J Twi 
F(z) = = —+ + O(z), 
XB. y2z Y By? 2 
toe se 
so that 


Py(t) = 2] 4 — 4 
v ‘ 
; Biy; 





No 
a 


ke ee ee Dw ee 


where Y(t) represents the contribution of the other 
poles of P,,(z). These poles are the non-vanishing 
roots of F(z) = 0. If we let z = s*, this equation 
can be written 
: lil, Ky(ajsIy(yjs) ~ I (ajs)Ky(yjs) 
g=1 yl Kol ajsy(y;s) + In(aj;s)Ky(y;s) 





0. (15) 


It can be shown that the roots of this equation are 
conjugate complex and lie on the imaginary s axis. 

Letting z, k = 1, 2, .... denote the non-zero 
roots of F(z) = 0 and letting /z = ix, z =x? so 
that the xz are roots of 





G(x) = 
Y1 (a;x) Jy (y;x) — Jy 0 a;x)Y, (y; x) 
& iis 1(4;x)], (7; J 1(ajx)¥4 (y; =0, (16) 
j=1 Yolajx)Jy (yjx) — Jol ajx) Vy (jx) 
we have 
, co otk? 
y(t)= 23 .(17) 


k=1 x G’(x,) 
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Evaluating this expression with the aid of the 
eigenvalue Eq. 16, we obtain 








en *kt 
Yitj=2- 23 
4 k=1 ‘ ; 1 Taff 43, ‘ + $4.) 
j=1 , $6, jk 
bb . (18) 
where 
Do, jk = Yol asp) IY XR) — In (ajxg) Yy (yj xp) 


1 jk = Yj (a5xp) J (yj; %,) - Jy (a; xp) ¥y (yj xp), 


This expression for Y(t) is correct for all values 
of t. In the calculations performed, it was found 
that the times for which the solution was evaluated 
were large enough that only the smaller roots xz 
were needed to evaluate Y(t). 

For such values of t, Eqs. 16 and 18 were sim- 


plified to 











n 14 (y;x) 
Gx)= & > By =0 (19) 
- = Ii(yjx) Ind ax — Y1(y;x) 
and 
Y(t) = 
ied 
.%§ me 
4 k=1 1 — J ?(y; xp) 
ha? - 
] Ph (yjxR) In aj, _ Yi (yj m) 
- (20) 


respectively. The exces ad the approximations 
checks in each case calculated. 

If we assume that the first root of Eq. 19 is 
small, further approximations can be made for the 
calculation of this root. For the case of two layers, 
Eq. 19 then becomes 


2 F X 
p 4. In ot Yj 


It 
Oo 





+ ypx 3 
Li a ee 


From this, the smallest root x, can a approximated 


" 2(y7By hi ¥2Bq) 


oe "mn. 3 i: 

7373 |Pa(in2- 2) «8, (in - 2)| 
(21) 
It is interesting to note that this expression 
coincides with the exponent obtained in the sim- 

plified theory of Tempelaar-Lietz for two layers. 
The transform of the production rate from the jth 

layer is given by 


x2 = 











W312) 
j% 
Gj(z) = ~~ : 
i z 2 hy (2) 
jat Bja; Wo jz) 
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By a procedure similar to the one used in obtaining 
P,,(t), we obtain 


aft) Bi; 








= Yi |) se Oo nn tee wee cP ar? 
q =Byp 7 ~e 
J] 
where 
Z,(t) = 
amy Ge” ae 
9 = ik 
“Fi at 1 x2 (92 + $F ) 
k= ~@? 
kt 
$ B; 4 0,7 1,jk 
=I é2 
u 0,7k 
. (23) 


The xz are the non-vanishing roots of Eq. 16. 
Again for large values of t, a simplified expres- 
sion for Eq. 23 can be constructed. This is 











. J (7%) Pos 
c w16%j*k) In 2 ajxp ~ Yy (yj xR) 
Tat 1 — J iy; Xp) 
py = Bia 2 
Gh ]1(¥j%») In 5 ax, - ¥, (yu) 
. (24) 


A auutinn treatment isis che: solution P; (r,t). 


THE PRESSURE AT SMALL TIMES 


The use of the solutions derived in the preceding 
section is not convenient for small times, since 
too many terms in the expressions for Y(t) and 
Z;(t) are required. However, since it is known that 
a bounded reservoir behaves as an infinite reservoir 
during its early life, the solutions corresponding 
to infinite and bounded reservoirs must coincide 
for all times smaller than the time at which the 
influence of the boundary is said to be felt upon 
the behavior of the reservoir. 

The actual procedure used in the calculation of 
the results obtained in this paper was, first, to 
evaluate the solution corresponding to an infinite 
field, and, then, to evaluate the solution corre- 
sponding to the bounded reservoir, using enough 
terms in the series for the transient parts of the 
solutions to obtain a significant time range where 
the two solutions would overlap. 


The solution to the problem of finding the well 
pressure in an infinite field composed of n horizon- 
tal layers has been solved by Horner.* The equa- 
tions to be solved are the same as for the bounded 
reservoir, but the boundary equation in Eq. 3 must 
be replaced by the condition 


P(r,t) +0 asr+o, t>0, 
hat 2 Tee Seen aa 


The solution for P,,(t) can be written 





*The authors gratefully acknowledge this contribution which 
has only been published within the Shell companies. 



























An 
728; oP 
where * eae _ and 
ne 4B 
J(t) =f :}. as 
u m2 
n B; 
+ y2 y2 
j= fa: : 
0j 0j : it. 
Yn; Y. . 
(Se lorleig ty)" 6 3 oe 
=1 Jo; +9; mu \j=1 eit; 
+ (27) 


The following abbreviations have been used. 
Joj = Jo(ajVu) , 
Jag = In(ajVu) , 
Yoj = Yo(a; Vu) ’ 
Yj; = ¥; (a;Vu) . 

The integral in the form given in Eq. 27 is quite 
intractable for numerical calculations. For the case 
of n = 2 (that is, when the reservoir is composed 
of only two layers), however, approximations can 
be obtained which make a numerical analysis 
possible. 

Introducing approximations for the Bessel func- 
tions for small values of u, we can show that 





~ I(t) = ByBy1n*(a3/a3) 0 et du 
(B, +B)? 0 ul(Inu-A)* +27] 
Terre Terr ST . (28) 
where 
_ Yr 
By In“; a3 + Bp In y at 





B, + Bo 


Now if we let r = te4, the integral in Eq. 28 is 
transformed into the integral 


e~ur d 
dr) = f - 


u(In? u+72) 
This integral must be evaluated to calculate J(t). 
From the table of Laplace transforms of McLach- 
lan, Humbert and Poli> it can be found that 





oo udu 
pest. fo . 
adele. OS? 


This last integral can be further transformed as 
follows. 


rt” du 


co 7 dy co omtl 
es T(u+1)~ 


alla+)D” co 


oo x ut+nd | 
> pes sets Fi 
n=0 0 [(ui+n+1) H inte 





ntu ‘ 


F(r,u) = > —_——_—_———__ 
1 n=0 IT'(u+n+1) 


This function F(z,u) satisfied the differential 
equation 
dF lia pu-l 
dr I(u) 
with the condition F(0,u) = 0. Solving this equation 
we find that 





F(7,u) = “UV yu-l dy = r(i- tee) 
7, U Wau tat ev’y y=e Tu) 
where 
co 
Our) = fe v¥-! dy 


T 
is an incomplete gamma function. Then, upon sub- 
stitution of these results, 
(u,7) 
(r) =e? pd du. 
d f T(u ) . . . . . . . . (29) 
This expression is used to calculate (7) for 


values of r< 10. For values of 7 > 10, an asymptotic 
expression is used for Q(u,7) which gives rise to 
the equation 





1 dv l-ri v dv 
— J r*D'(1-v) +? o TMI (1 -v) . 
2 
4/ ae. . Ge) 


25 r*T(1 -—v) 
It is found that the nenien of the last two terms 
decreases rapidly with increasing values of r. 

An upper limit for ¢{7) can be obtained in the 
following manner. 


l dv © dv iL 
< —< — =-e 
(7) J 7 J 7 


Inr 


This upper limit gives a good approximation for 
(7) for large values of r. 
The values of G(r) obtained are given in Table 1. 





TABLE 1 ~ VALUES OF ¢(7) 








; gir) : Plt) 
1 0.452 1x 104 0. 1003 
2 0.392 2x 104 0.0939 
3 0.359 4x 104 0.0883 
4 0.337 1x 105 0.0818 
6 0.308 2x 105 0.0775 

1x 10 0.2765 4x 10° 0.0736 

2x 10 0.2398 1x 106 0.0690 

4x 10 0.2100 2x 10° 0.0659 

1x 102 0.1792 4x 10° 0.0631 

2x 10? 0. 1608 1x 107 0.0597 

4x 10" 0.1456 

1x 10 0.1291 

2x 10 0.1189 

4x 103 0.1101 
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THE PRODUCTION RATE AT SMALL TIMES 


The production rate at small times is obtained 
from the solution for an infinite reservoir evaluated 
in the case of two layers. The Laplace transform 
of the fraction production rate from the first layer 
is 








K,(a, V2) 
az) 4 11 Kofay v2) 
q 2 Ky(a, v2) Ky (ag yz) 
~~... 2 Pee 
Ko( a; Vz) Ko(a2 Vz) 
é bike Se Mie ae @ wm ae % . .G)) 


The corresponding solution is 
q(t) 1 BH? qy(z) 
qm Be 4 

the path of integration passing to the right of the 


origin. Replacing the Bessel functions by the ap- 
proximations 


e! dz , 





Koa; Vz) = - Int a,vz , 
a 1 
Ky(a;Vz) = 


? 
a, Vz 
and neglecting quantities of the order a,, it can be 
shown that 





q(t) 
‘ = 
1 B Fad By ln ¥ azvz et 
Qi im y a... 
B-i B, In a2 V2 + By In5 a, vz 
or i 
q(t) = 2 “7 By f 3 


= ———_—_| e**az, 


ee ee ee ee ee ee ‘ a 632) 
where A is defined as before. The first term in the 
integrand yields the constant 6; /(8; + 82). The 
integral 

1 B+ ez? dz 


sone —  e —— ee te ee (3 
2ni ale z(1n z —A) " 


is evaluated as follows. According to Doetsch,® 


aoe (34) 
z(In z-—A) 
is the Laplace transform of 
os) 7 A 
J Tus dus (r = te”) ~ (35) 


The function Eq. 34 has a pole at z = eA which, 
however, does not contribute to the integral Eq. 33. 
This contribution, therefore, must be subtracted, 
yielding 
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a B+i°° et dz co 7% dy . 
2ni ee z(In z — A) SPa+w ane 





We can then write 


[Bro 


a 
oe » rr, 
9 By+Bq (By +B)? 


It can be noted that values of ¢{r) have already 
been calculated for evaluation of the pressures. 

The effects of the transient terms in the infinite 
solutions are both proportional to ¢(7); but, whereas 
this effect is negligible in the pressure equation 
for all but very small times, the effect of this term 
on the rate solution is in general of considerable 
magnitude. 


mb 





d(r) . .(36) 


INCLUSION OF THE SKIN EFFECT 


The skin factor S; of the jth layer is defined by 
the equation * 

a 

i Qf t) p,/ 2mk sb ; . phd 

The equations to be solved are the same in the 
presence of a skin as without one, except that the 
boundary condition (that the pressures at the well 
are the same in all layers) must be restated to say 
that the pressures inside the wellbore are the same 
in all layers. This can be stated as 

S59j(t)p; 


P;.(t) = P,,(t) 


a ae 
Pee) 


where P,,(t) is independent of }. 

The method of analysis used to solve this prob- 
lem is identical to the method outlined for solving 
the equations without the skin; the equations are 
similar, and only some extra terms containing Sj 
are present. 

The transform of P,,(t) can be obtained as 


— 1 
P(e) 0 . . .(39) 
] wb; (z) 





where 


W(z) = aS; Vz yz) + Uofz) , 


and Wp (z) and W1;(z) have the same meaning as 
as before. Similarly 


Wy (2) — 


, q;(z) = 2nBja; Vz Oo  — 


vj (z) 


The functions P,,,(z) and 9j(z) are again meromor- 
phic functions of z, and the inversion proceeds in 
the same manner as before. Then, 








*Where Py; denotes the pressure drop at the well in the for- 
mation in L 


ayer j. 














+ Ys(t) , (41) 


where Y,;(t) represents the transient solution. 

With the same approximations as before and in 
the time range of validity of these approximations, 
the skin factors S; always appear in such a fashion 
that they can be completely absorbed by replacing 
a; by a’, which is defined by 


a 4 eS 


The solution for the infinite reservoir with skins 
was found by a method similar to the one used by 
Horner; it was concluded that again, in the time 
range of validity of the approximations taken when 
no skin was present, the presence of the skin 
factors in the equations could be absorbed into 
the quantities a; in the manner of Eq. 42. 

Hence,,. it can be concluded that, for times of 
practical importance, the effect of a skin with skin 
factor S; in each one of the layers can be described 
completely by assigning a fictitious well radius 
Owe to each layer, the relation between the true 
well radius r,,,; and toss being given by 


* 


-S; 
en oo. el | Re S)) 


NOMENCLATURE 


= recovery/unit pressure drop of Layer j 
(Eq. 14) 
thickness of jth layer, cm 
h1 + bo, total thickness, cm 
productivity index of Layer j (Eq. 15) 


permeability of jth layer, darcy 
(kyhy = kogho)/ (hy + ho ds weighted mean 
of permeabilities, darcy 
= pressure in the jth layer at position r 
and time ¢t, atm 
= pressure at the well at time ¢, atm 
= average pressure in the jth layer at 
time ¢ 
p* = pressure obtained by straight-line ex- 
trapolation of first linear portion of 


Pw vs In[AtXt + At)] 


P fj = pressure in the formation of Layer 7 at 
the well 


P;(r,t) = pressure drop in the jth layer at posi- 
tion 7 and time ¢ 


P,,(t) = pressure drop function at the well 


total production rate from multilayer 
reservoir in cc/sec, measured at ini- 
tial reservoir conditions 

total production rate in cc/sec measured 
at average reservoir conditions 

production rate from jth layer at time ¢, 
in cc/sec, measured at initial reser- 
voir conditions 

flow rate out of jth layer at shut-in time 
At for a well closed in at production 
time t 

production rate during after-production 
at dimensionless shut-in time At 

cumulative production from Layer j at 
time t, in cc, measured at initial res- 
ervoir conditions 

fictitious well radius to represent skin 
effect, cm 

skin factor of jth layer, dimensionless 


production time, days 

time at which the pressure disturbance 
reaches the point 7, sec 

transient solution in well-pressure equa- 
tion 

transient solution in rate equation for 
jth layer 

loading constant (dimensionless) with 
respect to T 

differential depletion per unit thickness 
at time t between the layers of a two- 
layer reservoir 

shut-in time, sec 

dimensionless shut-in time = kAt/dycrZ 

oil-filled porosity of the jth layer, di- 
mensionless 

(hyhy + hoh2)/(h1 + h2) weighted mean 
of porosities 

Euler’s constant, In y = 0.5772 
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Complete Your Library of Permanent Literature With 


Petroleum Transactions Covering Years 1925 through 1960 


CONTENT OF VOLUMES ORDER FORM 
Reprinting of AIME Petroleum 
Transactions has been accomplished INDICATE NUMBER OF BOOKS DESIRED IN COLUMN SHOWN 
by combining two of the original vol- PRICE 
umes into one book for Volumes G-25 hg nna ——— 
through 179. Material on production Original Price Number Price Number | Add 
statistics and certain chapters of non- Fee Pr wl i fee tak bork. tate 
technical material were omitted in the Saal 
reprinting. The books in this series, , . 6.30 
containing two of the vohanes, 1925-1926 G-25—G-26 Book $9.00 $6. 
are six by nine inches or 1927-28-29  77- 82 $9.00 $6.30 
Volumes 186 through 198 were re- 1930-31 86- 92 $9.00 $6.30 
printed in their original form, with one 1932.33 98-103 $9.00 $6.30 
book of 812 by 10% inches in size 
containing 2 single volume. 1934-35 107-114 $9.00 $6.30 
1936-37 118-123 $9.00 $6.30 


BINDING, PACKAGING, PRICES 1938-39 127-132 $9.00 $6.30 


Books are bound in standard AIME 1940-41 136-142 $9.00 $6.30 
red cloth binding, and individually 1942-43 146-151 $9.00 $6.30 
packaged in cardboard cartons for de- : 
livery. Prices are listed on the order _—- 155-160 $9.00 $6.30 
form at right for cash purchase of 1946-47 165-170 $9.00 $6.30 
eg books, — a mre of 1948-49 174-179 $9.00 $6.30 

0 per cent to E mem as 
shown. No provision is made for 1949 186 Large Book $7.00 $4.90 
credit or installment purchases, or for 1950 189 $7.00 $4.90 
discount on purchase of complete set 1951 192 $7.00 $4.90 
of books. D ill be postpaid. 
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1953 198 $7.00 $4.90 
BONUS BOOKS 
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the Institute through 1952. Prices are 1980 219 $7.00 $4.90 
$5.00 to Non-Members and $3.50 to Index (1921-1953) $5.00 $3.50 
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Analyzer 


‘ 


~ Stream flow lines shown in white, signal paths in black 


” 


iis block diagram shows how six process 
reams can be selectively analyzed by the 
E Model 184-B. For brevity, Stream #1 
ily will be followed through the system. 
Stream #1 is tapped, and the sample 
vaporized in the Liquid Sample Vaporizer 
\). Sample #1 then flows (with the other 
ive vaporized samples) into the enclosed, 
ieated Multi-Stream Solenoid Switching 
Manifold (G). The solenoids are electric- 
ally directed by the Programmer (D). As 
it is selected, Sample #1 is valved into the 


Final Sampling Assembly (B), where it is 
filtered and metered into the Analyzer. The 
Analyzer quantitatively and qualitatively 
separates the sample’s components by gas 
chromatography, produces signals which 
appear as a bar graph on Recorder (E). 

The signal can also be fed into a Pneu- 
matic Control Accessory (F) which pro- 
duces a 15 psi signal for automatic control 
of one variable component of the stream. 

The Analyzer with Reverse Flow Valve 
(C) allows for the analysis of light com- 


ponents individually plus a “total heavies 
measurement; it back-flushes and analyzes 
heavy ends as a single peak, shortening 
analysis time in the process. The Program- 
mer (D) is available in a variety of modes 
for multi-component and multi-stream 
analysis. It controls the Switching Mani- 
fold (G) and the Analyzer, and the read- 
out of the Recorder (E). 

_All flow lines and sampling building 
blocks are kept above the dew point, where 
necessary, by steam tracing. 


UNTIL NOW—PROCESS CONTROL INSTRUMENTATION 
LIKE THIS HAD TO BE CUSTOM-BUILT 


Standard building blocks give the new P-E Model 184-B Process Vapor Fractometer the 


4 flexibility of the most advanced special engineering developments in process gas chromatography. 


A new, highly versatile instrument, the. Model 184-B 
brings automatic process control by gas chromatog- 
raphy within the reach of virtually any petroleum or 
chemical processing plant. 

The Model 184-B uses building block sampling and 
Sensing units in various combinations to produce proc- 
ess analytical systems which mean continuing extra 
dis idends through reduced maintenance, low-cost train- 
inv and highest reliability. 

‘hese standard building blocks include multi-port 


sampling valves for backflushing and analyzing heavy 
ends ... versatile multi-column systems ... multi- 
stream programmer kits... pneumatic units for auto- 
matic control .... and a wide variety of sampling devices. 

Perkin-Elmer, the world’s largest manufacturer of 
gas chromatographic instruments, will help you in your 
process instrumentation requirements. For more infor- 
mation, write for a new descriptive booklet. 


All building blocks in any P-E system. conform to standard 
refinery and chemical plant safety practices. 


INSTRUMENT DIVISION 


Perkin-Elmer Gyno 


NORWALK, 


CONNECTICUT 
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GIVES: YOU DOWN-TO-EARTH RESULT 
SPACE—the Bendix G-20 Automatic Programming ‘‘package’’— sets new standards fbr 
ease of use, power and efficiency. Designed in concert with the G-20 computer syste 

SPACE is complemented by numerous advanced equipment features. Here are the auto- 
matic programming methods which form an integral part of SPACE... and Bendix G-20 


systems, large or small: 





SPAR —Symbolic Assembly Programming. Allows the programmer to maintain direct 
Introduction to control over all G-20 operations. Provides the efficiency of machine language program. 
Programming Systema | : : Pek: 
; ming without the complexities. 





AL.COM —an Algebraic Compiler based on the international notation of ALGOL. Easy: 
to-use ALCOM permits the statement of scientific problems in natural mathematical 
language...simplifies and speeds problem solving. 


COBOL _—Common Business Oriented Language permits statement of data processing 
problems in natural business language for high-speed computer solution... makes flex: 
i} e){-UKX- We) mr-1) ©) py-] ol-1 4 (ommmel-Loll net-] bur- [ale Ms) o\-\el[-] Imei at-1¢-[ea-1 acm 


EXECUTIVE — provides automatic program scheduling and component assignment 
...permits maximum-efficiency in parallel processing and utilization of components 


See for yourself how SPACE...combined with outstanding equipment capabilities ... 
has put the G-20 in a class by itself. Investigate today. For your copy of ‘‘Introduction‘to 
G-20 Programming Systems,” write, wire or call: 


Condy” Bendix Computer Division 











